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Abstract 

The understanding of fluid turbulence has considerably progressed in recent years. 
The application of the methods of statistical mechanics to the description of the 
motion of fluid particles, i.e. to the Lagrangian dynamics, has led to a new quantita- 
tive theory of intcrmittency in turbulent transport. The first analytical description 
of anomalous scaling laws in turbulence has been obtained. The underlying physi- 
cal mechanism reveals the role of statistical integrals of motion in non-equilibrium 
systems. For turbulent transport, the statistical conservation laws are hidden in 
the evolution of groups of fluid particles and arise from the competition between 
the expansion of a group and the change of its geometry. By breaking the scale- 
invariance symmetry, the statistically conserved quantities lead to the observed 
anomalous scaling of transported fields. Lagrangian methods also shed new light 
on some practical issues, such as mixing and turbulent magnetic dynamo. 
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"Well," said Pooh, " we keep looking for Home and not finding it, so 1 thought 
that if we looked for this Pit, we'd be sure not to find it, which would be a Good 
Thing, because then we might find something that we weren't looking for, which 
might be just what we were looking for, really" . A. Milne, Tigger is unbounced. 
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I. INTRODUCTION 



Turbulence is the last great unsolved problem of classical physics which has evaded 
physical understanding and systematic description for many decades. Turbulence is a state 
of a physical system with many degrees of freedom strongly deviating from equilibrium. The 
first obstacle to its understanding stems from the large number of degrees of freedom actively 
involved in the problem. The scale of injection, where turbulence is excited, usually differs 
dramatically from the scale of damping, where dissipation takes place. Nonlinear interactions 
strongly couple the degrees of freedom by transferring excitations from the injection to the 
damping scale throughout the so-called inertial range of scales. The ensuing complicated 
and irregular dynamics calls for a statistical description. The main physical problem is to 
understand to what extent the statistics in the inertial interval is universal, i.e. independent 
of the conditions of excitation and dissipation. In such general formulation, the issue goes 
far beyond fluid mechanics, even though the main examples and experimental data are 
provided by turbulence in continuous media. From the standpoint of theoretical physics, 
turbulence is a non-equilibrium field-theoretical problem with many strongly interacting 
degrees of freedom. The second deeply rooted obstacle to its understanding is that far 
from equilibrium we do not possess any general guiding rule, like the Gibbs principle in 
equilibrium statistical physics. Indeed, to describe the single-time statistics of equilibrium 
systems, the only thing we need is the knowledge of dynamic integrals of motion. Then, our 
probability distribution in phase space is uniform over the surfaces of constant integrals of 
motion. Dynamically conserved quantities play an important role in turbulence description, 
too, as they flow throughout the inertial range in a cascade-like process. However, the 
conserved quantity alone does not allow one to describe the whole statistics but only a 
single correlation function which corresponds to its flux. The major problem is to obtain 
the rest of the statistics. 

In every case, the starting point is to identify the dynamical integral of motion that 
cascades through the inertial interval. Let us consider the forced 3d Navier-Stokes equation 
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d t v{r, t) + v(r, t) ■ V»(r, f) - uV 2 v(r, t) = -Vp{r, t) + f(r, t), 



(1) 



supplemented by the incompressibility condition V • v — 0. An example of injection mecha- 
nism is a random large-scale forcing f(r, t) with correlation length L. The relevant integral of 
motion, conserved in the absence of injection and dissipation, is kinetic energy / v 2 dr/2 and 
the quantity which cascades throughout the inertial interval is energy density in wavenumber 
space. The energy flux- constancy relation was derived in Kolmogorov (1941) and it involves 
the third-order moment of the longitudinal velocity increments: 



The separation r is supposed to lie in the inertial interval, ranging from the injection scale 
L down to the viscous dissipation scale. The major physical assumption made to derive the 
so-called 4/5 law is that the mean energy dissipation rate e v = u((Vv) 2 ) has a nonzero limit 
as the viscosity v tends to zero. This clearly points to the non-equilibrium flux nature of 
turbulence. The assumption of finite dissipation gives probably the first example of what is 
called "anomaly" in modern field-theoretical language: A symmetry of the inviscid equation 
(here, time-reversal invariance) is broken by the presence of the viscous term, even though 
the latter might have been expected to become negligible in the limit of vanishing viscosity. 
Note that the 4/5 law (2) implies that the third-order moment is universal, that is, it depends 
on the injection and the dissipation only via the mean energy injection rate, coinciding with 
e v in the stationary state. To obtain the rest of the statistics, a natural first step made by 
Kolmogorov himself was to assume the statistics in the inertial range be scale invariant. The 
scale invariance amounts to assuming that the probability distribution function (PDF) of the 
rescaled velocity differences r~ h A r v can be made r-independent for an appropriate h. The 
n-th order moment of the longitudinal velocity increments ((A r t>) n ) (structure functions) 
would then depend on the separation as a power law r CTn with the "normal scaling" behavior 
a n = hn. The rescaling exponent may be determined by the flux law, e.g. h — 1/3 for 3d 
Navier-Stokes turbulence. In the original Kolmogorov theory, the scale invariance was in fact 
following from the postulate of complete universality: the dependence on the injection and 
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the dissipation is carried entirely by e v not only for the third-order moment but for the whole 
statistics of the velocity increments. The velocity difference PDF could then involve only 
the dimensionless combination (e 1 ,r)~ 1//3 A r t> and would be scale invariant. There are cases, 
like weakly nonlinear wave turbulence (Zakharov et al, 1992), where both scale-invariance 
and complete universality are assured by the fact that the statistics in the inertial range 
is close to Gaussian. That does not hold for strongly nonlinear systems. Already in 1942, 
L. D. Landau pointed out that all the velocity structure functions (except the third one) are 
averages of nonlinear functions of the flux. They are therefore sensitive to its fluctuations, 
which depend on the specific injection mechanisms. Consequently, the velocity statistics in 
the inertial range may have nonuniversal features. 

Experiments do not support scale invariance either. The structure functions are in fact 
found experimentally to have a power-law dependence on the separation r. However, the 
PDF of the velocity differences at various separations cannot be collapsed one onto another 
by simple rescaling and the scaling exponent a n of the structure functions is a nonlinear 
concave function of the order n. As the separation decreases in the inertial range, the 
PDF becomes more and more non-Gaussian, with a sharpening central peak and a tail 
that becomes longer and longer. In other words, the smaller the separations considered, 
the higher the probability of very weak and strong fluctuations. This manifests itself as a 
sequence of strong fluctuations alternating with quiescent periods, which is indeed observed 
in turbulence signals and is known as the phenomenon of intermittency. The violation of 
the dimensional predictions for the scaling laws is referred to as "anomalous scaling" for it 
reflects, again, a symmetry breaking. The Euler equation is scale-invariant and the scales 
of injection and dissipation are supposed to be very large and small (formally, the limits to 
infinity and zero should be taken). However, the dynamics of turbulence is such that the 
limits are singular and scale invariance is broken. The presence of a finite injection scale L, 
irrespective of its large value, is felt throughout the inertial range precisely via the anomalies 
{(A r v) n ) oc (e v r) n / 3 (L/r) n/3_<T ". 

The non-Gaussianity of the statistics, the anomalous scaling and the intermittency of 
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the field occur as a rule rather than exception in the context of fluid dynamics. The same 
phenomenology is observed in many other physical systems. An incomplete list includes 
compressible Navier-Stokes turbulence, Burgers' turbulence, scalar and magnetic fields. Ex- 
amples of scalar fields are provided by the temperature of a fluid, the humidity in the 
atmosphere, the concentration of chemical or biological species. The advection-diffusion 
equation governing the transport of a nonreacting scalar field by an incompressible velocity 
is: 



where <p describes the sources. For scalar dynamics, the space integral of any function of 
9 is conserved in the absence of sources and diffusion. In their presence, the corresponding 
relation for the flux of 9 2 was derived in Yaglom (1949): 



The major physical assumption is again that the mean scalar dissipation rate e = k((V9) ) 
remains finite even in the limit where the molecular diffusivity k vanishes. Consider the par- 
ticular case when the advecting velocity v satisfies the 3d Navier-Stokes equation. Assuming 
again scale-invariance, the flux relations (2) and (4) would imply that the scaling exponent 
of both the velocity and the scalar field is 1/3. As it was expected for the velocity, the 
scalar structure functions S n (r) = ([0(r, t) — 9(0, t)] n ) would then depend on the separation 
as power laws r^ n with ( n = n/3. Experiments indicate that scale invariance is violated for 
a scalar field as well, that is ( n ^ n/3. More importantly, the intermittency of the scalar is 
much stronger than that of the velocity, in particular, n/3 — ( n is substantially larger than 
n/3 — cr n . It was a major intuition of R.H. Kraichnan to realize that the passive scalar could 
then be intermittent even in the absence of any intermittency of the advecting velocity. 

The main ambition of the modern theory of turbulence is to explain the physical mech- 
anisms of intermittency and anomalous scaling in different physical systems, and to under- 
stand what is really universal in the inertial-interval statistics. It is quite clear that strongly 



d t 9{r, t) + v(r, t) • V0(r, f) - kV 2 9(t, t) = ip(r, t), 



(3) 
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non-equilibrium systems generally do not enjoy the same degree of universality as those in 
equilibrium. In the absence of a unified approach to non-equilibrium situations, one tries 
to solve problems on a case-by-case basis, with the hope to learn if any universal guiding 
principle may be recognized. It is in solving the particular problems of passive scalar and 
magnetic fields that an important step in general understanding of turbulence has been re- 
cently made. The language most suitable for the description of the systems turned out to 
be the Lagrangian statistical formalism, i.e. the description of the motion of fluid particles. 
This line of analysis, pioneered by L.F. Richardson and G.I. Taylor in the twenties and 
later developed by R.H. Kraichnan and others, has been particularly effective here. The 
results differ from case to case. Some fields are non- Gaussian but scale invariance is not 
broken, while others have turned out to be amenable to the first ever analytical description 
of anomalous scaling laws. The anomalous exponents have been found to be universal, but 
not the constants appearing in the prefactors of generic correlation functions. This has 
provided a quantitative clarification of Landau's previously mentioned remark and of the 
aspects of turbulence statistics that may still be expected to be universal. More importantly, 
the anomalous scaling has been traced to the existence of statistical integrals of motion. The 
mechanism is quite robust and relevant for transport by generic turbulent flows. The na- 
ture of those integrals of motion strongly differs from that of the dynamic conservation laws 
that determine equilibrium statistics. For any finite number of fluid particles, the conserved 
quantities are functions of the interparticle separations that are statistically preserved as the 
particles are transported by the random flow. For example, at scales where the velocity field 
is spatially smooth, the average distance R between two particles generally grows exponen- 
tially, while the ensemble average (R~ d ) is asymptotically time-independent in a statistically 
isotropic <i-dimensional random flow. The integrals of motion change with the number of 
particles and generally depend nontrivially on the geometry of their configurations. In the 
connection between the advected fields and the particles, the order of the correlation func- 
tions is equal to the number of particles and the anomalous scaling issue may be recast as 
a problem in statistical geometry. The nonlinear behavior of the scaling exponents with the 
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order is then due to the dependence of the integrals of motion on the number of particles. 
The existence of statistical conservation laws signals that the Lagrangian dynamics keeps 
trace of the particle initial configuration throughout the evolution. This memory is what 
makes the correlation functions at any small scale sensitive to the presence of a finite injec- 
tion length L. We believe that, more generally, the notion of statistical integrals of motion 
is a key to understand the universal part of the steady-state statistics for systems far from 
equilibrium. 

The aim of this review is a description of fluid turbulence from the Lagrangian viewpoint. 
Classical literature on Lagrangian dynamics mostly concentrated on turbulent diffusion and 
pair dispersion, i.e. the distance traveled by one particle or the separation between two 
particles as a function of time. By contrast, in that general picture that has emerged recently, 
the evolution of the multiparticle-configuration geometry takes center stage. The main body 
of the review will present these novel aspects of Lagrangian dynamics and their consequences 
for the advected fields. We shall adhere to the following plan. The knowledge accumulated 
on one and two particle dynamics has been extensively covered in literature (Pope, 1994; 
Majda and Kramer, 1999). The objective of the first three Sections in Chapter II is to point 
out a few fundamental issues, with particular attention to the basic differences between the 
cases of spatially smooth and nonsmooth velocity fields. We then proceed to the multiparticle 
statistics and the analysis of hidden statistical conservation laws that cause the breakdown 
of scale-invariance. Most of this analysis is carried out under the assumption of a prescribed 
statistics of the velocity field. In Chapter III we shall analyze passive scalar and vector 
fields transported by turbulent flow and what can be inferred about their statistics from 
the motion of fluid particles. In Chapter IV, we briefly discuss the Lagrangian dynamics 
in the Burgers and the Navier-Stokes equations. The statistics of the advecting velocity 
is not prescribed anymore, but it results from nonlinear dynamics. Conclusions focus on 
the impact of the results presented in this review on major directions of future research. 
Readers from other fields of physics interested mainly in the breakdown of scale invariance 
and statistical conservation laws may restrict themselves to Sects. II.C, II.E, III.C, V. 
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The picture presented in this review is, to a large extent, an outcome of joint work and 
numerous discussions with our colleagues, E. Balkovsky, D. Bernard, A. Celani, M. Chertkov, 
G. Eyink, A. Fouxon, U. Frisch, I. Kolokolov, A. Kupiainen, V. Lebedev, A. Mazzino and 
A. Noullez. We thank K. Khanin, P. Muratore-Ginanneschi, A. Shafarenko, B. Shraiman and 
the referee for valuable comments about the manuscript. We are indebted to R. H. Kraichnan 
whose works and personality have been a permanent source of inspiration. 

II. PARTICLES IN FLUID TURBULENCE 

As explained in the Introduction, understanding the properties of transported fields 
involves the analysis of the behavior of fluid particles. We have therefore decided to first 
present results on the time-dependent statistics of the Lagrangian trajectories R n (t) and to 
devote the subsequent Chapter III to the description of transported fields. In the present 
Chapter we sequentially increase the number of particles involved in the problem. We start 
from a single trajectory whose effective motion is a simple diffusion at times longer than 
the velocity correlation time in the Lagrangian frame (Sect. II. A). We then move to two 
particles. The separation law of two close trajectories depends on the scaling properties of 
the velocity field v{r,t). If the velocity is smooth, that is \v(R n ) — v(R m )\ oc \R n — R m \, 
then the initial separation grows exponentially in time (Sect. II. B). The smooth case can 
be analyzed in much detail using the large deviation arguments presented in Sect. II.B.l. 
The reader mainly interested in applications to transported fields might wish to take the 
final results (21) and (27) for granted, skipping their derivation and the analysis of the 
few solvable cases where the large deviations may be calculated exactly. If the velocity 
is nonsmooth, that is \v(R n ) — v(R m )\ oc \Rn — Rm\ a with a < 1, then the separation 
distance between two trajectories grows as a power of time (Sect. II. C), as first observed 
by Richardson (1926). We discuss important implications of such a behavior on the nature 
of the Lagrangian dynamics. The difference between the incompressible flows, where the 
trajectories generally separate, and compressible ones, where they may cluster, is discussed 
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in Sect. II. D. Finally, in the consideration of three or more trajectories, the new issue of 
geometry appears. Statistical conservation laws come to light in two-particle problem and 
then feature prominently in the consideration of multiparticle configurations. Geometry and 
statistical conservation laws are the main subject of Sect. II. E. Although we try to keep the 
discussion as general as possible, much of the insight into the trajectory dynamics is obtained 
by studying simple random ensembles of synthetic velocities where exact calculations are 
possible. The latter serve to illustrate the general features of the particle dynamics. 

A. Single-particle diffusion 

The Lagrangian trajectory R(t) of a fluid particle advected by a prescribed incompress- 
ible velocity field v(r,t) in d space dimensions and undergoing molecular diffusion with 
diffusivity k is governed by the stochastic equation (Taylor, 1921), customarily written for 
differentials: 



Here, (3(t) is the ci-dimensional standard Brownian motion with zero average and covariance 
function {(3 l (t) f3' J '(f)) = 5 %] min(£, f). The solution of (5) is fixed by prescribing the particle 
position at a fixed time, e.g. the initial position R(0). 

The simplest instance of (5) is the Brownian motion, where the advection is ab- 
sent. The probability density V(AR; t) of the displacement AR(t) = R(t) — R(0) sat- 
isfies the heat equation (d t — kV 2 )V = whose solution is the Gaussian distribution 
V(AR; t) = (AnKf)~ d l 2 exp[— (AR) 2 /(Ant)]. The other limiting case is pure advection with- 
out noise. The properties of the displacement depend then on the specific trajectory under 
consideration. We shall always work in the frame of reference with no mean flow. We assume 
statistical homogeneity of the Eulerian velocities which implies that the Lagrangian velocity 
V(t) = v(R(t),i) is statistically independent of the initial position. If, additionally, the 




(5) 
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Eulerian velocity is statistically stationary, then so is the Lagrangian one 1 . The single-time 
expectations of the Lagrangian velocity coincide in particular with those of the Eulerian 
one, e.g. (V(t)) = (v) = 0. The relation between the multi-time statistics of the Eulerian 
and the Lagrangian velocities is however quite involved in the general case. 
For k — 0, the mean square displacement satisfies the differential equation: 



d_ 

dt 



(AR(t)) 2 ) = 2 /' (V(t) ■V(s))ds = 2f t (V(0) • V(s)) ds , (6) 



where the second equality uses the stationarity of V(t). The behavior of the displacement 
is crucially dependent on the range of temporal correlations of the Lagrangian velocity. Let 
us define the Lagrangian correlation time as 

= S~(V(0).V(8))da 

(V 2 ) ' 1 1 

The value of r provides a measure of the Lagrangian velocity memory. Divergence of r 
is symptomatic of persistent correlations. As we shall discuss in the sequel, no general 
relation between the Eulerian and the Lagrangian correlation times can be established but 
for the case of short-correlated velocities. For times f « r, the 2-point function in (6) 
is approximately equal to (v 2 ) and the particle transport is ballistic: ((AR) 2 } ~ (v 2 )t 2 . 
When the Lagrangian correlation time is finite, a generic situation in a turbulent flow, 
an effective diffusive regime arises for t > t with ((AR) 2 ) = 2(v 2 } rt (Taylor, 1921). The 
particle displacements over time segments spaced by distances much larger than r are indeed 
almost independent. At long times, the displacement AR behaves then as a sum of many 
independent variables and falls into the class of stationary processes governed by the Central 
Limit Theorem. In other words, the displacement for t ^> r becomes a Brownian motion in 
d dimensions with 



1 This follows by averaging the expectations involving V(t + r) over the initial position R(0) 
(on which they do not depend) and by the change of variables R(0) > R(t) under the velocity 
ensemble average. The argument requires the incompressibility of the velocity, see Sect. II. D. 
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(AiT(t) AR?(t)) ~ 2Z#'t, (8) 

where 

oo 

= i J (V\0) V j (s) + V j (0) V l (s)} ds . (9) 
o 

The same arguments carry over to the case of a non-vanishing molecular diffusivity. The 
symmetric second order tensor D l J describes the effective diffusivity (also called eddy diffu- 
sivity). The trace of D l J is equal to the long-time value (v 2 )t of the integral in (6), while its 
tensorial properties reflect the rotational symmetries of the velocity field. If it is isotropic, 
the tensor reduces to a diagonal form characterized by a single scalar value. The main 
problem of turbulent diffusion is to obtain the effective diffusivity tensor, given the velocity 
field v and the value of the diffusivity k. Exhaustive reviews of the problem are available 
in the literature (Bensoussan et ai, 1978; Pope, 1994; Fannjiang and Papanicolaou, 1996; 
Majda and Kramer, 1999). 

The other general issue in turbulent diffusion is about the conditions on the velocity 
field ensuring the Lagrangian correlation time r be finite and an effective diffusion regime 
take place for large enough times. A sufficient condition (Kraichnan, 1970; Avellaneda and 
Majda, 1989; Avellaneda and Vergassola, 1995) is that the vector potential variance (A 2 ) is 
finite, where the 3d incompressible velocity v = V x A. Similar conditions are valid for any 
space dimension. The condition k ^ is essential to the validity of the previous result, as 
shown by the counter-example of Rayleigh-Benard convective cells, see e.g. (Normand et ai, 
1977). In the absence of molecular noise, the particle circulates forever in the same convective 
cell, with no diffusion taking place at any time. This provides an example of subdiffusion : 
the integral in (6) goes to zero as t — > oo and the growth of the mean square displacement 
is slower than linear. Note that any finite molecular diffusivity, however small, creates thin 
diffusive layers at the boundaries of the cells; particles can then jump from one cell to 
another and diffuse. Subdiffusion is particularly relevant for static 2c? flows, where tools 
borrowed from percolation/statistical topography find most fruitful applications (Isichenko, 
1992). Trapping effects required for subdiffusion are, generally speaking, favored by the 
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compressibility of the velocity field, e.g., in random potentials (Bouchaud and Georges, 
1990). Subdiffusive effects are expected to be overwhelmed by chaotic mixing in flows 
leading to Lagrangian chaos, i.e., to particle trajectories that are chaotic in the absence of 
molecular diffusion (Ottino, 1989; Bohr et al, 1998). This is the generic situation for 3d 
and 2d time-dependent incompressible flows. 

Physical situations having an infinite Lagrangian correlation time r correspond to su- 
perdiffusive transport : divergences of the integral in (6) as t — > oo signal that the particle 
transport is faster than diffusive. A classical example of such behavior is the class of par- 
allel flows presented by Matheron and de Marsily (1980). If the large-scale components of 
the velocity field are sufficiently strong to make the particle move in the same direction for 
arbitrarily long periods the resulting mean square displacement grows more rapidly than t. 
Other simple examples of super diffusive motion are Levy-type models (Geisel et al, 1985; 
Shlesinger et al, 1987). A detailed review of super diffusive processes in Hamiltonian systems 
and symplectic maps can be found in Shlesinger et al. (1993). 

Having listed different subdiffusive and super diffusive cases, from now on we shall be 
interested in random turbulent flows with finite Lagrangian correlation times, which are 
experimentally known to occur for sufficiently high Reynolds numbers (Pope, 1994). For the 
long-time description of the diffusion in such flows, it is useful to consider the extreme case 
of random homogeneous and stationary Eulerian velocities with a short correlation time. 
The formal way to get these processes is to change the time scale by taking the scaling limit 
lim /i^v(r,/it), i.e. considering the process as viewed in a sped-up film. We assume that 
the connected correlation functions 2 decay fast enough when time differences increase. The 
elementary consequences of those assumptions are the existence of the long-time asymptotic 
limit and the fact that it is governed by the Central Limit Theorem. When fi — > oo, we 

2 The connected correlation functions, also called cumulants, are recursively defined by the relation 
{vi ... v n ) = J2 II ((^(l), • • • > v Tv a (n a ))) w ith the sum over the partitions of {1, . . . ,n}. 



13 



recover a velocity field which is Gaussian and white in time, characterized by the 2-point 
function 



The advection by such velocity fields was first considered by Kraichnan (1968) and it is 
common to call the Gaussian ensemble of velocities with 2-point function (10) the Kraichnan 
ensemble. For the Kraichnan ensemble, the Lagrangian velocity V(t) has the same white 
noise temporal statistics as the Eulerian one v(r,t) for fixed r and the displacement along 
a Lagrangian trajectory AR(t) is a Brownian motion for all times. The eddy diffusivity 
tensor is D l i = D lj (0), which is a special case of relation (9). In the presence of molecular 
diffusion, the overall diffusivity is the sum of the eddy contribution and the molecular value 



In realistic turbulent flows, the Lagrangian correlation time r is comparable to the char- 
acteristic time scale of large eddies. Progress in numerical simulations (Yeung, 1997) and 
experimental technique (Voth et al, 1998; La Porta et al, 2001; Mordant et al. 2001) has 
provided information on the single particle statistics in the regime intermediate between bal- 
listic and diffusive. Such behavior is captured by the the subtracted Lagrangian autocorrela- 
tion function (V(0)(V(0) — V(t))) or its second time derivative that is the autocorrelation 
function of the Lagrangian acceleration. This information has provided stringent tests on 
simple stochastic models (that eliminate velocity fields), often used in the past to describe 
the one-particle and two-particle statistics in turbulent flows (Pope, 1994). The Kraichnan 
ensemble that models stochastic velocity fields, certainly missrepresents the single particle 
statistics by suppressing the regime of times smaller than r. It constitutes, however, as we 
shall see in the sequel, an important theoretical laboratory for studying the multiparticle 
statistics in fluid turbulence. 




(10) 



K5 ij . 
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B. Two-particle dispersion in a spatially smooth velocity 

The separation i? 12 — R\ — R 2 between two fluid particles with trajectories R n {t) = 
R(t; r n ) passing at t — through the points r n satisfies (in the absence of Brownian motion) 
the equation 

R 12 = v(R 1 ,t)-v(R 2 ,t). (11) 

We consider first an incompressible flow where the particles generally separate. In this 
Section, we start from the smallest distances where the velocity field can be considered 
spatially smooth due to viscous effects. In next Section, we treat the dispersion problem 
for larger distances (in the inertial interval of turbulence) where the velocity field has a 
nontrivial scaling. Finally, we describe a compressible flow and show how the separation 
among the particles is replaced by their clustering as the degree of compressibility grows. 

1. General considerations 

In smooth velocities, for separations R 12 much smaller than the viscous scale of tur- 
bulence, i.e. in the so-called Batchelor regime (Batchelor, 1959), we may approximate 
v(R 1 ,t) — v(R 2 ,t) a(t) Rn(t) with the Lagrangian strain matrix <T^(t) = V ' jV l (R 2 (i),t). 
In this regime, the separation obeys the ordinary differential equation 

R 12 (t) = a(t) R 12 (t) , (12) 

leading to the linear propagation 

Rn(t) = W(t)R 12 (0), (13) 

where the evolution matrix is defined as W l ^(t) = dR l (r\i) / dr^ with r = r 2 . We shall also 
use the notation W(t; r) when we wish to keep track of the initial point or W(t; r, s) if the 
initial time s is different from zero. 

The equation (12), with the strain treated as given, may be explicitly solved for arbitrary 
a(t) only in the Id case by expressing W(t) as the exponential of the time-integrated strain: 
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ln[R(t)/R(0)} = In W(t) = j a{s) ds = X 



Jo 



(14) 



We have omitted subscripts replacing R u by -R. When t is much larger than the correlation 
time r of the strain, the variable X behaves as a sum of many independent equally distributed 
random numbers X = Y^i Vi with iV oc t/r. Its mean value (X) = N(y) grows linearly in 
time. Its fluctuations X — (X) on the scale 0{t x l 2 ) are governed by the Central Limit 
Theorem that states that (X — (X))/N l l 2 becomes for large N a Gaussian random variable 
with variance (y 2 ) — (y) 2 = A. Finally, its fluctuations on the larger scale 0(t) are governed 
by the Large Deviation Theorem that states that the PDF of X has asymptotically the form 



This is an easy consequence of the exponential dependence on N of the generating function 
(e zX ) of the moments of X. Indeed, (e zX ) = e NS ( z \ where we have denoted (e zy ) = e s ^ 
(assuming that the expectation exists for all complex z). The PDF V(X) is then given by 
the inverse Laplace transform ^ / e~ zX+NS ( z ) dz with the integral over any axis parallel to 
the imaginary one. For X oc N, the integral is dominated by the saddle point z Q such that 
S'(zo) = X/N and the large deviation relation (15) follows with H = —S(z ) + z S'(z ). 
The function H of the variable X/N — (y) is called entropy function as it appears also in the 
thermodynamic limit in statistical physics (Ellis, 1985). A few important properties of H 
(also called rate or Cramer function) may be established independently of the distribution 
V(y). It is a convex function which takes its minimum at zero, i.e. for X taking its mean 
value (X) = NS'(0). The minimal value of H vanishes since S(0) = 0. The entropy is 
quadratic around its minimum with H"(0) = A -1 , where A = S"(0) is the variance of y. 
The possible non-Gaussianity of the y } s leads to a non-quadratic behavior of H for (large) 
deviations of X/N from the mean of the order of A/S""(0). 

Coming back to the logarithm \nW(t) of the interparticle distance ratio in (14), its 
growth (or decay) rate A = (X)/t is called the Lyapunov exponent. The moments ([R(t)] n ) 
behave exponentially as exp[7(n)t] with 7(71) a convex function of n vanishing at the origin. 



V(X) oc e 



-NH(X/N-{y)) 



(15) 
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Even if A = 7'(0) < 0, high-order moments of R may grow exponentially in time, see, for 
instance, the behavior of the interparticle distance discussed in Section II. D. In this case, 
there must be one more zero ri\ of 7(71) and a statistical integral of motion, (R ni ), that does 
not depend on time at large times. 

In the multidimensional case, the solution (13) for R(t) is determined by products of 
random matrices rather than just random numbers. The evolution matrix W(t) may be 
written as 



W(t) = Texp 



1 t OO 3 2 

Ja(s)ds = J ' cr(s n )ds n . . . J a(s 2 ) ds 2 J cr(si) ctei . (16) 
n=0 00 

This time-ordered exponential form is, of course, not very useful for direct calculations except 

for the particular case of a short-correlated strain, see below. The main statistical properties 

of the separation vector R needed for most physical applications might still be established 

for quite arbitrary strains with finite temporal correlations. The basic idea goes back to 

Lyapunov (1907) and Furstenberg and Kesten (1960) and it found further development in 

the Multiplicative Ergodic Theorem of Oseledec (1968). The modulus R of the separation 

vector may be expressed via the positive symmetric matrix W T W. The main result states 

that in almost every realization of the strain, the matrix |lniy T iy stabilizes as t — > 00. 

In particular, its eigenvectors tend to d fixed ortho normal eigenvectors f \. To understand 

that intuitively, consider some fluid volume, say a sphere, which evolves into an elongated 

ellipsoid at later times. As time increases, the ellipsoid is more and more elongated and 

it is less and less likely that the hierarchy of the ellipsoid axes will change. The limiting 

eigenvalues 

\= limrHnlWfil (17) 

t^oo 

define the so-called Lyapunov exponents. The major property of the Lyapunov exponents 
is that they are realization-independent if the strain is ergodic. The usual convention is to 
arrange the exponents in non-increasing order. 

The relation (17) tells that two fluid particles separated initially by -R(O) pointing into 
the direction f i will separate (or converge) asymptotically as exp(Ajt). The incompressibility 
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constraints det(VF) = 1 and — imply that a positive Lyapunov exponent will exist 
whenever at least one of the exponents is nonzero. Consider indeed 

E(n) = lim r 1 In (\R(t) / R(0)} n ) , (18) 

whose slope at the origin gives the largest Lyapunov exponent Ai. The function E{n) ob- 
viously vanishes at the origin. Furthermore, E(—d) = 0, i.e. incompressibility and isotropy 
make that (R~ d ) is time-independent as t — > oo (Furstenberg, 1963; Zeldovich et al., 1984). 
Negative moments of orders n < — 1 are indeed dominated by the contribution of directions 
R(0) almost aligned to the eigendirections f 2 , ■ ■ ■ fa- At n < 1 — d the main contribution 
comes from a small subset of directions in a solid angle oc exp(dA^) around f d . It follows 
immediately that (R n ) oc exp[\ d (d + n)t\ and that (R~ d ) is a statistical integral of motion. 
Since E{n) is a convex function, it cannot have other zeroes except — d and if it does not 
vanish identically between those values. It follows that the slope at the origin, and thus 
Ai, is positive. The simplest way to appreciate intuitively the existence of a positive Lya- 
punov exponent is to consider, following Zel'dovich et al. (1984), the saddle-point 2d flow 
v x = Xx,v y = —\y. A vector initially forming an angle <f> with the x-axis will be stretched 
after time T if cos0 > [1 + exp(2AT)]~ 1//2 , i.e. the fraction of stretched directions is larger 
than 1/2. 

A major consequence of the existence of a positive Lyapunov exponent for any random 
incompressible flow is an exponential growth of the interparticle distance R(t). In a smooth 
flow, it is also possible to analyze the statistics of the set of vectors R(t) and to establish 
a multidimensional analog of (15) for the general case of a nondegenerate Lyapunov ex- 
ponent spectrum. The final results will be the Large Deviation expressions (21) and (27) 
below. The idea is to reduce the <i-dimensional problem to a set of d scalar problems ex- 
cluding the angular degrees of freedom. We describe this procedure following Balkovsky and 
Fouxon (1999). Consider the matrix I(t) = W(t)W T (t) ) representing the tensor of inertia 
of a fluid element like the above mentioned ellipsoid. The matrix is obtained by averaging 
R l (t)R j (t)/£ 2 d over the initial vectors of length I. In contrast to W T W that stabilizes at 
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large times, the matrix / rotates in every realization. To account for that rotation, we rep- 
resent the matrix as T AO with the orthogonal O composed of the eigenvectors of I and 
the diagonal A having the eigenvalues e 2pi , . . . e 2f>d arranged in non-increasing order. The 
evolution equation d t I — al + Ia T takes then the form 



d t pi = a ii) a = OaO T , (19) 

a,o = ao, n ii = e -^±y^, (20) 

with no summation over repeated indices. We assume isotropy so that at large times the 
SO(d) rotation matrix O is distributed uniformly over the sphere. Our task is to describe 
the statistics of the stretching and the contraction, governed by the eigenvalues pi. We 
see from (19,20) that the evolution of the eigenvalues is generally entangled to that of the 
angular degrees of freedom. As time increases, however, the eigenvalues will become widely 
separated (pi ... p d ) for a majority of the realizations and fijj — > a^ for i < j 
(the upper triangular part of the matrix follows from antisymmetry). The dynamics of 
the angular degrees of freedom becomes then independent of the eigenvalues and the set of 
equations (19) reduces to a scalar form. The solution pi = / ' au(s) ds allows the application 
of the large deviation theory, giving the asymptotic PDF: 

V(p!,...,p d ;t) oc exp[-tH(p 1 /t-\ 1 ,...,p d ^ 1 /t-\ d ^ 1 )] 

x 9( Pl - p 2 ) . . . Q(p d _ x - p d ) 5( Pl + . . . + p d ) . (21) 

The Lyapunov exponents \ are related to the strain statistics as \ = (an) where the 
average is temporal. The expression (21) is not valid near the boundaries pi = p i+ \ in a 
region of order unity, negligible with respect to \{t at times t ^> (Aj — Aj + i) _1 . 

The entropy function H depends on the details of the strain statistics and has the same 
general properties as above: it is non-negative, convex and it vanishes at zero. Near the 
minimum, H(x) m ^{C^ v )ijXiXj with the coefficients of the quadratic form given by the 
integrals of the connected correlation functions of a defined in (19): 

dj = J ((vn(t),Vjj(t'))) dt' , i,j = l,...,d-l. (22) 
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In the 5-correlated case, the entropy is everywhere quadratic. For a generic initial vector 
r, the long-time asymptotics of ln(i?/r) coincides with V(pi) = J V(pi, . . . , pa) dp 2 . . . dpd 
which also takes the large- deviation form at large times, as follows from (21). The quadratic 
expansion of the entropy near its minimum corresponds to the lognormal distribution for 
the distance between two particles 



with r = R(0), A = Ai and A = C n . 

It is interesting to note that under the same assumption of nondegenerate Lyapunov 
spectrum one can analyze the eigenvectors of the evolution matrix W (Goldhirsch et ai, 
1987). Note the distinction between the eigenvectors of W and f i of W T W. Let us 
order the eigenvectors according to their eigenvalues. Those are real due to the assumed 
nondegeneracy and they behave asymptotically as exp(Ait), . . . , exp(A^t). The eigenvector 
converges exponentially to a fixed vector and any subspace spanned by {e d „ fc , . . . , e d } for 
< k < d tends asymptotically to a fixed subspace for every realization. Remark that the 
subspace is fixed in time but changes with the realization. 

Molecular diffusion is incorporated into the above picture by replacing the differential 
equation (12) by its noisy version 



The separation vector is subject to the independent noises of two particles, hence the factor 
2 with respect to (5). The solution to the inhomogeneous linear stochastic equation (24) 
is easy to express via the matrix W(t) in (16). The tensor of inertia of a fluid element 
P j (t) = ^R l (t)R j (t) is now averaged both over the initial vectors of length £ and the 
noise, thus obtaining (Balkovsky and Fouxon, 1999): 

t 

I(t) = W{t)W{t) T + ^ d J W(t) W{8) T W{8)\- x W{t) T d8. (25) 



The matrix I(t) evolves according to d t I = a I + la + ^ and the elimination of the 




(23) 



dR(t) = a(t)R(t) dt + 2y/Rdp{t) . 



(24) 
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angular degrees of freedom proceeds as previously. An additional diffusive term 2k exp(— 2pi) 
appears in (19) and the asymptotic solution becomes 

pi(t) = au(s)ds + ^ In |l + ^ J exp - 2 J a u (s') ds' dsj . (26) 

The last term in (26) is essential for the directions corresponding to negative \. The 
molecular noise will indeed start to affect the motion of the marked fluid volume when the 
respective dimension gets sufficiently small. If t is the initial size, the required condition 
Pi ~ —pi = — \n(i 2 \\i\/ k) is typically met for times t ~ p*/\\i\. For longer times, the 
respective pi is prevented by diffusion to decrease much below — p*, while the negative Aj 
prevents it from increasing. As a result, the corresponding pi becomes a stationary random 
process with a mean of the order — p*. The relaxation times to the stationary distribution 
are determined by a, which is diffusion independent, and they are thus much smaller than t. 
On the other hand, the components pj corresponding to non-negative Lyapunov exponents 
are the integrals over the whole evolution time t. Their values at time t are not sensitive to 
the latest period of evolution lasting over the relaxation time for the contracting pi. Fixing 
the values of pj at times t ^> p*/\\i\ will not affect the distribution of the contracting pi and 
the whole PDF is factorized (Shraiman and Siggia, 1994; Chertkov et al, 1997; Balkovsky 
and Fouxon, 1999). For Lagrangian dynamics in 3d developed Navier-Stokes turbulence 
there are, for instance, two positive and one negative Lyapunov exponents (Girimaji and 
Pope, 1990). For times t > P3/I A 3 | we have then 

V oc exp[-tH( Pl /t-\ 1 ,p 2 /t-\ 2 )} V st (p 3 ), (27) 

with the same function H as in (21) since P3 is independent of p\ and p 2 . Note that the 
account of the molecular noise violates the condition Pi — as fluid elements at scales 



smaller than ^Jk,/\\ 3 \ cannot be distinguished. To avoid misunderstanding, note that (27) 
does not mean that the fluid is getting compressible: the simple statement is that if one 
tries to follow any marked volume, the molecular diffusion makes this volume statistically 
growing. 
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Note that we have implicitly assumed i to be smaller than the viscous length r\ = y^/|A 3 | 
but larger than the diffusion scale ^Jk/\X 3 \. Even though v and k are both due to molecular 
motion, their ratio widely varies depending on the type of material. The theory of this 
section is applicable for materials having large Schmidt (or Prandtl) numbers v/k. 

The universal forms (21) and (27) for the two-particle dispersion are basically everything 
we need for physical applications. We will show in the next Chapter that the highest 
Lyapunov exponent determines the small-scale statistics of a passively advected scalar in a 
smooth incompressible flow. For other problems, the whole spectrum of exponents and even 
the form of the entropy function are relevant. 

2. Solvable cases 

The Lyapunov spectrum and the entropy function can be derived exactly from the given 
statistics of a for few limiting cases only. The case of a short-correlated strain allows for a 
complete solution. For a finite-correlated strain, one can express analytically A and A for 
a 2d long-correlated strain and at large space dimensionality. 

i) Short-correlated strain. Consider the case where the strain a(t) is a stationary 
white-in-time Gaussian process with zero mean and the 2-point function 

K(f)<7 M (0) = 28(t-lf)C ijkt . (28) 

This case may be considered as the long-time scaling limit lim ^a(fit) of a general strain 
along a Lagrangian trajectory, provided its temporal correlations decay fast enough. It may 
be also viewed as describing the strain in the Kraichnan ensemble of velocities decorrelated 
in time and smooth in space. In the latter case, the matrix C^i = —VjVeD tk (0), where 
D lj (r) is the spatial part in the 2-point velocity correlation (10). We assume D lj (r) to 
be smooth in r (or at least twice differentiate), a property assured by a fast decay of its 
Fourier transform D l i{k). Incompressibility, isotropy and parity invariance impose the form 
£>«(r) = D 5 ij -\d' l i(r) with 
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tfJ(r) = L>i[(d+l)^'r 2 



2r v] + o(r 2 ). 



(29) 



The corresponding expression for the 2-point function of a reads 



Cijkt = D 1 [(d + l)SikSje — 5ij5ke — SuSjk] , 



(30) 



with the constant Di having the dimension of the inverse of time. 

The solution of the stochastic differential equation (12) is given by (16) with the matrix 
W(t) involving stochastic integrals over time. For a white-correlated strain, such integrals are 
not defined unambiguously but require a regularization that reflects finer details of the strain 
correlations wiped out in the scaling limit. An elementary discussion of this issue may be 
found in the Appendix. For an incompressible strain, however, the ambiguity in the integrals 
defining W(t) disappears so that we do not need to care about such subtleties. The random 
evolution matrices W(t) form a diffusion process on the group SL(d) of real matrices with 
unit determinant. Its generator is a second-order differential operator identified by Shraiman 
and Siggia (1995) as M = D\[dH 2 — (d+ 1)</ 2 ], where H 2 and J 2 are the quadratic Casimir 
of SL(d) and its SO(d) subgroup. In other words, the PDF of W(t) satisfies the evolution 
equation (d t — M)V(W; t) = 0. The matrix W(t) may be viewed as a continuous product of 
independent random matrices. Such products in continuous or discrete versions have been 
extensively studied (Furstenberg and Kesten, 1960; Furstenberg, 1963; Le Page, 1982) and 
occur in many physical problems, e.g. in Id localization (Lifshitz et ai, 1988; Crisanti et ai, 



If we are interested in the statistics of stretching-contraction variables only, then 
W(t) may be projected onto the diagonal matrix A with positive non-increasing entries 
e 2pl , . . . ,e 2pd by the decomposition W = OA f O', where the matrices O and O' belong to 
the group SO(d). As observed in (Bernard et ai, 1998; Balkovsky and Fouxon, 1999), the 
generator of the resulting diffusion of pi is the rf-dimensional integrable Calogero-Sutherland 
Hamiltonian. The pi obey the stochastic Langevin equation 



1993). 



d t pi = Ad^coth(pi - pj) + r], 



(31) 
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where rj is a white noise with 2-point function (r]i(t) r)j(t')) = 2Di(d5ij — 1) S(t — t'). At long 
times the separation between the pj's becomes large and we may approximate coth(p^) by 
±1. It is then easy to solve (31) and find the explicit expression of the PDF (21): 

H ( x ) = ih E x i > A * = D < d - 2< + 1) ■ (32) 
i=i 

Note the quadratic form of the entropy, implying that the distribution of R{t) takes the 
lognormal form (23) with A = Ai and A = 2Di(d — 1). The calculation of the long-time 
distribution of the leading stretching rate p\ goes back to Kraichnan (1974). The whole 
set of d Lyapunov exponents was first computed in (Le Jan, 1985), see also (Baxendale, 
1986). Gamba and Kolokolov (1996) obtained the long-time asymptotics of pi by a path 
integral calculation. The spectral decomposition of the Calogero-Sutherland Hamiltonian, 
see (Olshanetsky et ai, 1983), permits to write explicitly the PDF of pi for all times, 
ii) 2d slow strain. In 2d, one can reduce the vector equation (12) to a second-order scalar 
form. Let us indeed consider the case of a slow strain satisfying & <C a 2 and differentiate 
the equation R = crR with respect to time. The term with & is negligible with respect to 
a 2 and a little miracle happens here: because of incompressibility, the matrix a is traceless 
and a 2 is proportional to the unit matrix in 2d. We thus come to a scalar equation for the 
wave function ^ = R x + iR y : 

d 2 V = (a 2 u + a 12 a 21 )q. (33) 

This is the stationary Schrodinger equation for a particle in the random potential U = 
S 2 — Q 2 , where S 2 = o\ x + (cr 12 + 02i) 2 /4 and Q 2 = {a 12 — cr 2 i) 2 /4 is the vorticity. Time 
plays here the role of the coordinate. Our problem is thus equivalent to localization in the 
quasi-classical limit (Lifshitz et ai, 1988) and finding the behavior of (33) with given initial 
conditions is similar to the computation of a Id sample resistivity, see e.g. (Abrikosov and 
Ryzhkin, 1978; Kolokolov, 1993). Based on these results we can assert that the modulus 
l^l = R in random potentials grows exponentially in time, with the same exponent that 
controls the decay of the localized wave function. 
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The problem can be solved using semi-classical methods. The flow is partitioned in 
elliptic (Q > S) and hyperbolic (S > Q) regions (Weiss, 1991), corresponding to classical 
allowed (U < 0) and forbidden (U > 0) regions. The wave function \1/ is given by two 
oscillating exponentials or one decreasing and one increasing, respectively. Furthermore, the 
typical length of the regions is the correlation time r, much larger than the inverse of the 
mean strain and vorticity S^ s and Q^ms- ^ follows that the increasing exponentials in the 
forbidden regions are large and dominate the growth of R(t). With exponential accuracy 
we have: 



where the real part restricts the integration to the hyperbolic regions. The parameters A 
and A in the lognormal expression (23) are immediately read from (34): 



Note that the vorticity gave no contribution in the ^-correlated case. For a finite correlation 
time, it suppresses the stretching by rotating fluid elements with respect to the axes of 
expansion. The real part in (35) is indeed filtering out the elliptic regions. Note that the 
Lyapunov exponent is given by a single-time average, while in the 5-correlated case it was 
expressed by the time-integral of a correlation function. It follows that A does not depend on 
the correlation time r and it can be estimated as S rms for Q rms ^ S rms . The corresponding 
estimate of the variance is A ~ (S 2 )t . As the vorticity increases, the rotation takes over, 
the stretching is suppressed and A reduces. One may show that the correlation time r s of 
the stretching rate is the minimum between 1/Q rms and r (Chertkov et al, 1995a). For 
fl rms T > 1 we are back to a 5-correlated case and A ~ (S 2 ) /Q rms . All those estimates can 
be made systematic for a Gaussian strain (Chertkov et al, 1995a). 

iii) Large space dimensionality. The key remark for this case is that scalar products like 
R l (ti)R l (t2) are sums of a large number of random terms. The fluctuations of such sums 
are vanishing in the large-ci limit and they obey closed equations that can be effectively 




(34) 





(35) 
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studied for arbitrary strain statistics. This approach, developed in (Falkovich et al, 1998), 
is inspired by the large-iV methods in quantum field theory ('t Hooft, 1974) and statistical 
mechanics (Stanley, 1968). Here, we shall relate the behavior of the interparticle distance to 
the strain statistics and find explicitly A and A in (23). The strain is taken Gaussian with 
zero mean and correlation function {<Tij(t)<Jke(0)) = (2-D /rd) Sn-Sjg g(() where higher-order 
terms in 1/d are neglected. The integral of g is normalized to unity and ( = t/r. At large 
d, the correlation function F = (i?*(ti)i? 4 (t 2 )) satisfies the equation 



with the initial condition d^F((, 0) = 0. The limit of large d is crucial for the factorization 
of the average leading to the second equality in (36). The dimensionless parameter (3 = 2Dt 
measures whether the strain is long or short-correlated. Since (36) is linear and the coefficient 
on the right-hand side depends only on the time difference, the solution may be written as 
a sum of harmonics F x (Ci, C2) = exp[Ar(Ci + C2)] ^(Ci — C2)- Inserting it into (36), we get 
the Schrodinger equation for the even function ^(t): 



At large times, the dominant contribution comes from the largest exponent A corresponding 
to the ground state in the potential —@g(£)- Note that the "energy" is proportional to —A 2 
and A is the Lyapunov exponent since (R 2 ) = F ((,()■ From quantum mechanics textbooks 
it is known that the ground state energy in deep and shallow potentials is proportional to 
their depth and its square, respectively. We conclude that A oc D for a fast strain (small f3) 
and A oc \J D/t in the slow case (large (5). At large time differences, the potential term in 
(37) is negligible and A determines both the growth of (R 2 ) and the decay of different-time 
correlation function at Ci + C2 fixed. That also shows that the correlation function becomes 
independent of the larger of the times t\ and t 2 when their difference exceeds r. 

For the fast strain case, one can put g(() = 5(() and the solution of (37) is amazingly 
simple: F(d,(2) = R 2 (0) exp[/5min(Ci, (2)]- The Lyapunov exponent A = D, in agreement 



F(CiX2) = r 2 (a tJ (t 1 )a lk (t 2 )R\t 1 )R k (t 2 ))=f3g(Ci-C2)F(Ci,C2), 



(36) 



^(0 + [-(^) 2 + /5^(C)]^(C) = o. 



(37) 
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with the result Ai = Did 2 + 0(d) obtained for the Kraichnan ensemble. For the slow case, 
the stretching rate is independent of r at a given value of D/r (determining the simultaneous 
correlation function of the strain). The analysis of the Schrodinger equation (37) with a deep 
potential also gives the correlation time r s of the stretching rate, which does not generally 
coincide with the strain correlation time r (Falkovich et al, 1998). 

C. Two-particle dispersion in a nonsmooth incompressible flow 

In this Section we study the separation between two trajectories in the inertial range of 
scales rj <^ r <^ L. The scales n and L stand in a 3d turbulent flow for the viscous and the 
injection scales (the latter is also called integral scale in that case). For a 2d inverse energy 
cascade flow, they would stand for the scales of injection and friction damping respectively. 
We shall see that the behavior of the trajectories is quite different from that in smooth flows 
analyzed previously. 

1. Richardson law 

As discussed in the Introduction, velocity differences in the inertial interval exhibit 
an approximate scaling expressed by the power law behavior of the structure functions 
{(A r v) n ) oc r a ". Low-order exponents are close to the Kolmogorov prediction a n = an with 
a — 1/3. A linear dependence of a n on n would signal the scaling A r v oc r a with a sharp 
value of a. A nonlinear dependence of a n indicates the presence of a whole spectrum of 
exponents, depending on the space-time position in the flow (the so-called phenomenon of 
multiscaling). The 2d inverse and the 3d direct energy cascades in Navier-Stokes equation 
provide concrete examples of the two possible situations. Rewriting (11) for the fluid particle 
separation as R = Av(R,t), we infer that dR 2 /dt = 2 R ■ Av oc R 1+a . If the value of a 
is fixed and smaller than unity, this is solved (ignoring the space-time dependence in the 
proportionality constant) by 

R l ~ a {t) - ^""(O) oc t, (38) 
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implying that the dependence on the initial separation is quickly wiped out and that R 
grows as t 1 ^ 1- "). For the random process R(t), the relation (38) is, of course, of the mean 
field type and should pertain to the long-time behavior of the averages 

(R<(t)) oc t^- a K (39) 

That implies their superdiffusive growth, faster than the diffusive one oc t^ 2 . The scaling 
law (39) might be amplified to the rescaling property 

V(R;t) = XV(XR; X^t) (40) 

of the interparticle distance PDF. Possible deviations from a linear behavior in the order ( 
of the exponents in (39) should be interpreted as a signal of multiscaling of the Lagrangian 
velocity Av(R(t),t) = AV(t). 

The power-law growth (39) for ( = 2 and a = 1/3, i.e. (R(t) 2 ) oc t 3 , is a direct conse- 
quence of the celebrated Richardson dispersion law (Richardson, 1926), the first quantitative 
phenomenological observation in developed turbulence. It states that 

±(R 2 ) oc (RY 3 . (41) 

The law (41) seems to be confirmed also by later experimental data, see Chapter 24 of Monin 
and Yaglom (1979) and (Jullien et al, 1999), and by the numerical simulations (Zovari et 
al, 1994; Elliott Jr. and Majda, 1996; Fung and Vassilicos, 1998; Boffetta et al, 1998). 
The more general property of self-similarity (40) (with a = 1/3) has been observed in the 
inverse cascade of two-dimensional turbulence (Jullien et al, 1999; Boffetta and Sokolov, 
2000; Boffetta and Celani, 2000). It is likely that (41) is exact in that situation, while it 
may be only approximately correct in 3d, although the experimental data do not allow yet 
to test it with sufficient confidence. 

It is important to remark that, even assuming the validity of the Richardson law (41), it 
is impossible to establish general properties of the PDF V(R; t) such as those in Sect. II.B.l 
for the single particle PDF. The physical reason is easy to understand if one writes 
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^ = 2n((AV) 2 >, (42) 

similarly to (6) and (7). Here r t = f*(AV(t) • AV(s)) ds /((AV) 2 } is the correlation time 
of the Lagrangian velocity differences. If d(R 2 )/dt is proportional to (R 2 ) 2 ^ 3 and ((AV) 2 ) 
behaves like (R 2 ) 1 ^ 3 then r t grows as (R 2 ) 1 ^ 3 oc t, i.e. the random process AV(t) is 
correlated across its whole span. The absence of decorrelation explains why the Central 
Limit Theorem and the Large Deviation Theory cannot be applied. In general, there is 
no a priori reason to expect V(R] t) to be Gaussian with respect to a power of R either, 
although, as we shall see, this is what essentially happens in the Kraichnan ensemble. 

2. Breakdown of the Lagrangian flow 

It is instructive to contrast the exponential growth (18) of the distance between the tra- 
jectories within the viscous range with the power law growth (39) in the inertial range. In the 
viscous regime the closer two trajectories are initially, the longer time is needed to reach a 
given separation. As a result, infinitesimally close trajectories never separate and trajectories 
in a fixed realization of the velocity field are continuously labeled by the initial conditions. 
Small deviations of the initial point are magnified exponentially, though. This sensitive de- 
pendence is usually considered as the defining feature of the chaotic behavior. Conversely, 
in the inertial interval the trajectories separate in a finite time independently of their initial 
distance -R(O), provided the latter was also in the inertial interval. The speed of this sepa- 
ration may depend on the detailed structure of the turbulent velocities, including their fine 
geometry (Fung and Vassilicos, 1998), but the very fact of the explosive separation is related 
to the scaling behavior A r v oc r a with a < 1. For high Reynolds numbers the viscous scale 
i] is negligibly small, a fraction of a millimeter in the turbulent atmosphere. Setting it to 
zero (or equivalently the Reynolds number to infinity) is an appropriate abstraction if we 
want to concentrate on the behavior of the trajectories in the inertial range. In such a limit, 
the power law separation between the trajectories extends down to arbitrarily small dis- 
tances: infinitesimally close trajectories still separate in a finite time. This makes a marked 
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difference in comparison to the smooth chaotic regime, clearly showing that developed tur- 
bulence and chaos are fundamentally different phenomena. As stressed in (Bernard et ai, 

1998) , the explosive separation of the trajectories results in a breakdown of the determinis- 
tic Lagrangian flow in the limit Re — > oo, see also (Frisch et al, 1998; Gaw§dzki, 1998 and 

1999) . The effect is dramatic since the trajectories cannot be labeled anymore by the initial 
conditions. Note that the sheer existence of the Lagrangian trajectories R(t; r) depending 
continuously on the initial position r would imply that lim ( \R(t; r±) — R(t; r 2 )\^ ) = . 
That would contradict the persistence of a power law separation of the Richardson type for 
infinitesimally close trajectories. Remark also that the breakdown of the deterministic La- 
grangian flow does not violate the theorem about the uniqueness of solutions of the ordinary 
differential equation R = v(R,t). Indeed, the theorem requires the velocity to be Lipschitz 
in r, i.e. that A r v < 0(r). As first noticed by Onsager (1949), the velocities for Re = oo are 
actually only Holder continuous: A r v ~ 0(r a ) with the exponent a < 1 (in Kolmogorov's 
phenomenology a = 1/3). The simple equation x = \x\ a provides a classical example with 
two solutions x = [(1 — a)^ 1 ^ 1- ") and x — 0, both starting from zero, for the non-Lipschitz 
case a < 1. It is then natural to expect the existence of multiple Lagrangian trajectories 
starting or ending at the same point. Such a possibility was first noticed and exploited in a 
somewhat different context in the study of weak solutions of the Euler equations (Brenier, 
1989; Shnirelman, 1999). Does then the Lagrangian description of the fluid break down 
completely at Re = oo? 

Even though the deterministic Lagrangian description is inapplicable, the statistical de- 
scription of the trajectories is still possible. As we have seen above, probabilistic questions 
like those about the averaged powers of the distance between initially close trajectories 
still have well defined answers. We expect that for a typical velocity realization, one may 
maintain at Re = oo a probabilistic description of the Lagrangian trajectories. In partic- 
ular, objects such as the PDF p(r, s; R, 1 1 v) of the time t particle position R, given its 
time s position r, should continue to make sense. For a regular velocity with deterministic 
trajectories, 
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p(r,s;R,t\ v) = S(R - R(t;r, s)) , (43) 

where R(t; r, s) denotes the unique Lagrangian trajectory passing at time s through r. In 
the presence of a small molecular diffusion, equation (5) for the Lagrangian trajectories has 
always a Markov process solution in each fixed velocity realization, irrespective of whether 
the latter be Lipschitz or Holder continuous (Stroock and Varadhan, 1979). The resulting 
Markov process is characterized by the transition probabilities p(r, s; R,t\ v) satisfying the 
advection-diffusion equation 3 

(d t - V R ■ v(R, t) - kV 2 r ) p(r, s;R,t\ v) = 0, (44) 

for t > s. The mathematical difference between smooth and rough velocities is that in the lat- 
ter case the transition probabilities are weak rather than strong solutions. What happens if 
we turn off the molecular diffusion? If the velocity is Lipschitz in r, then the Markov process 
describing the noisy trajectories concentrates on the deterministic Lagrangian trajectories 
and the transition probabilities converge to (43). It has been conjectured in (Gaw§dzki, 
1999) that, for a generic Re = oo turbulent flow, the Markov process describing the noisy 
trajectories still tends to a limit when k — > 0, but the limit stays diffused, see Fig. 1. In other 
words, the transition probability converges to a weak solution of the advection equation 

(d t - V R ■ v(R, t)) p(r, s;R,t\v) = , (45) 

which does not concentrate on a single trajectory, as it was the case in (43). We shall then 
say that the limiting Markov process defines a stochastic Lagrangian flow. This way the 
roughness of the velocity would result in the stochasticity of the particle trajectories persist- 
ing even in the limit k — > 0. To avoid misunderstanding, let us stress again that, according 
to this claim, the Lagrangian trajectories behave stochastically already in a fixed realization 
of the velocity field and for negligible molecular diffusivities, i.e. the effect is not due to the 



3 For k > and smooth velocities, the equation results from the ltd formula generalizing (A5) 
applied to (43) and averaged over the noise. 
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molecular noise or to random fluctuations of the velocities. This spontaneous stochasticity 
of fluid particles seems to constitute an important aspect of developed turbulence. It is 
an unescapable consequence of the Richardson dispersion law and of the Kolmogorov-like 
scaling of velocity differences in the limit Re — > oo and it provides for a natural mechanism 
assuring the persistence of dissipation in the inviscid limit: lim z/(| Vv| 2 ) 7^ 0. 

v— »0 

3. The example of the Kraichnan ensemble 

The general conjecture about the existence of stochastic Lagrangian flows for generic 
turbulent velocities, e.g. for weak solutions of the incompressible Euler equations locally 
dissipating energy, as discussed by Duchon and Robert (2000), has not been mathematically 
proven. The conjecture is known, however, to be true for the Kraichnan ensemble (10), as 
we are going to discuss in this Section. 

We should model the spatial part of the 2-point function (10) so that it has proper 
scalings in the viscous and inertial intervals. This can be conveniently achieved by taking 
its Fourier transform 



where cf J (r) scales as in the inertial interval i] <C r <C L, as r 2 in the viscous range 
r <tC 77 and tends to 2D 5 tj at very large scales r » I. As we discussed in Sect. II. A, 
Do gives the single-particle effective diffusivity. Notice that -Do — ®(L^) indicating that 
turbulent diffusion is controlled by the velocity fluctuations at large scales of order L. On 
the other hand, dij(r) describes the statistics of the velocity differences and it picks up 
contributions of all scales. In the limits rj — > and L — > 00, it takes the scaling form: 




(46) 



with < £ < 2. In physical space, 



D^(r) = D ft - \&(r), 



(47) 



lim tf'-(r) = D ir t[(d-1 + Z)ft - £— ], 



(48) 
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where the normalization constant D\ has the dimensionality of length 2 ^ time^ 1 . 

For < £ < 2 and rj > 0, the typical velocities are smooth in space with the scaling 
behavior visible only for scales much larger than the viscous cutoff 77. When the cutoff 
is set to zero, however, the velocity becomes nonsmooth. The Kraichnan ensemble is then 
supported on velocities that are Holder-continuous with the exponent £/2 — 0. That mimics 
the major property of turbulent velocities at the infinite Reynolds number. The limiting 
case i = 2 describes the Batchelor regime of the Kraichnan model: the velocity gradients are 
constant and the velocity differences are linear in space. This is the regime that the analysis 
of Sect. II.B.2(i) pertains to. In the other limiting case £ = 0, the typical velocities are 
very rough in space (distributional). For any £, the Kraichnan velocities have even rougher 
behavior in time. We may expect that the temporal roughness does not modify largely the 
qualitative picture of the trajectory behavior as it is the regularity of velocities in space, 
and not in time, that is crucial for the uniqueness of the trajectories (see, however, below). 

For time-decorrelated velocities, both terms on the right hand side of the Lagrangian 
equation (5) should be treated according to the rules of stochastic differential calculus. 
The choice of the regularization is irrelevant here even for compressible velocities, see Ap- 
pendix. The existence and the properties of solutions of such stochastic differential equations 
were extensively studied in the mathematical literature for velocities smooth in space, see 
e.g. (Kunita, 1990). Those results apply to our case as long as 77 > both for positive or 
vanishing diffusivity. The advection-diffusion equation (44) for the transition probabilities 
also becomes a stochastic equation for white-in-time velocities. The choice of the convention, 
however, is important here even for incompressible velocities: the equation should be inter- 
preted with the Stratonovich convention, see Appendix. The equivalent Ito form contains 
an extra second-order term that amounts to the replacement of the molecular diffusivity 
by the effective diffusivity (D + k) in (44). The Ito form of the equation explicitly ex- 
hibits the contribution of the eddy diffusivity, hidden in the convention for the Stratonovich 
form. As pointed out by Le Jan and Raimond (1998 and 1999), the regularizing effect of 
D permits to solve the equation by iteration also for the nonsmooth case giving rise to 
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transition probabilities p(r, s;R,t\ v) defined for almost all velocities of the Kraichnan en- 
semble. Moreover, the vanishing diffusivity limit of the transition probabilities exist, defining 
a stochastic Lagrangian flow. 

The velocity averages over the Kraichnan ensemble of the transition probabilities 
p(r, s; R,t\ v) are exactly calculable. We shall use a formal functional integral approach 
(Chertkov, 1997; Bernard et al, 1998). In the phase space path integral representation of 
the solution of (44), 

t 

f T J[»P(T)-(r(T)-»(r(T),T)) + Kp 2 (T)]dT 

p(r,s;R,t\v) = / e « VpVr , (49) 

r(s)=r 
r(t)=R 

for Syt, the Gaussian average over the velocities is easy to perform. It replaces the ex- 

t 

ponent in (49) by =F J[ip(r) ■ r(r) + (D + k)p 2 (t)]cIt and results in the path integral 
represention of the heat kernel of the Laplacian for which we shall use the operator notation 
e |t-s|(_D +K)v 2 ^ y ,. j n th er words, the average of (49) is the solution of the heat equation 
(with diffusivity D + k) equal to d(R — r) at time s. The above calculation confirms then 
the result discussed at the end of Sect. II. A about the all-time diffusive behavior of a single 
fluid particle in the Kraichnan ensemble. 

In order to study the two-particle dispersion, one should examine the joint PDF of the 
equal-time values of two fluid particles averaged over the velocities 

( p(n, s; Ri, 1 1 v) p{r 2 , s;R 2 ,t\ v)) = V 2 (r 1 ,r 2 ; Ri, R 2 ;t - s) . (50) 

The latter is given for the Kraichnan ensemble by the heat kernel e'* _s '- M2 (r , 1 , r 2 ; Ri, R 2 ) 
of the elliptic second-order differential operator 

2 2 

m 2 = J2 ^'K-v)v,v r , +«E^„. (51) 

n,n'=l n=l 

In other words, the PDF V 2 satisfies the equation (d t —M 2 )V 2 = S(t—s)S(Ri—ri)S(R 2 —r 2 ) : 
a result which goes back to the original work of Kraichnan (1968). Indeed, the Gaussian 
expectation (50) is again easily computable in view of the fact that the velocity enters 
through the exponential function in (49). The result is the path integral expression 
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, T/( E[ip„(r)-r„(r) + K p2( T )]+ £ £>« (r„(r)-r„, (r)) p ni (r) p n , 3 (r)) dr 

e s n=1 l[Vp n Vr n (52) 

»'n(s)=i'n ™ 
r n (t)=H n 

for the heat kernel of M. 2 . 

Let us concentrate on the relative separation R = Ri — R 2 of two fluid particles at 
time t, given their separation r at time zero. The relevant PDF V(r; R; t) is obtained by 
averaging over the simultaneous translations of the final (or initial) positions of the particles. 
Explicitly, it is given by the heat kernel of the operator M. = (cf J (r)+2/t S 1 *) V r »V r j equal to 
the restriction of M 2 to the translationally invariant sector. Note that the eddy diffusivity 
D , dominated by the integral scale, drops out from M.. The above result shows that 
the relative motion of two fluid particles is an effective diffusion with a distance-dependent 
diffusivity tensor scaling like in the inertial range. This is a precise realization of the 
scenario for turbulent diffusion put up by Richardson (1926). 

Similarly, the PDF of the distance R between two particles is given by the heat kernel 
e^ M (r;R), where M is the restriction of M 2 to the homogeneous and isotropic sector. 
Explicitly, 

1 



(d - l)D 1 r d - 1+ t + 2Kr d ~ l d r (53) 



in the scaling regime and its heat kernel may be readily analyzed. In the Batchelor regime 
£ = 2 and for k — > 0, the heat kernel of M reproduces the lognormal distribution (23) with 
A = 2£>i(d - 1) and A = D l( i(d - 1), see Sect. II.B.2(i). 

The simple criterion allowing to decide whether the Markov process stays diffused as 
k — ► is to control the limit r — > of the PDF V(r;R;t) (Bernard et al, 1998). For 
smooth velocities, it follows from (23) that 

lim V(r;R;t) = 5(R). (54) 

K— >0 

In simple words, when the initial points converge, so do the endpoints of the process. Con- 
versely, for < £ < 2 we have 
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in the scaling limit r] = 0, L = oo. That confirms the diffused character of the limiting 
process describing the Lagrangian trajectories in fixed non-Lipschitz velocities: the endpoints 
of the process stay at finite distance even if the initial points converge. If we set the viscous 
cutoff to zero keeping L finite, the behavior (55) crosses over for R 3> L to a simple 
diffusion with diffusivity 2D : at such large distances the particle velocities are essentially 
independent and the single particle behavior is recovered. 

The stretched-exponential PDF (55) has the scaling form (40) for a = £ — 1 and implies 
the power law growth (39) of the averaged powers of the distance between trajectories. The 
PDF is Gaussian in the rough case £ = 0. Note that the Richardson law (R 2 (t)} oc t 3 is 
reproduced for £ = 4/3 and not for £ = 2/3 (where the velocity has the spatial Holder 
exponent 1/3). The reason is that the velocity temporal decorrelation cannot be ignored 
and the mean-field relation (38) should be replaced by R l -^ 2 {t) - i2 1_ ^ 2 (0) oc f3(t) with 
the Brownian motion (3(t). Since (3(t) behaves as t 1 / 2 , the replacement changes the power 
and indeed reproduces the large-time PDF (55) up to a geometric power-law prefactor. 
In general, the time dependence of the velocities plays a role in determining whether the 
breakdown of deterministic Lagrangian flow occurs or not. Indeed, the relation (42) implies 
that the scale-dependence of the correlation time r t of the Lagrangian velocity differences 
may change the time behavior of (-R 2 ). In particular, (R 2 ) ceases to grow in time if r t oc (R 2 ) 13 
and ((AV) 2 ) oc (R 2 ) a with (3 > 1 — a. It has been recently shown in (Fannjiang et al, 2000) 
that the Lagrangian trajectories are deterministic in a Gaussian ensemble of velocities with 
Holder continuity in space and such fast time decorrelation on short scales. The Kolmogorov 
values of the exponents a — (5 — 1/3 satisfy, however, (5 < 1 — a. 

Note the special case of the average (R 2 ~^~ d ) in the Kraichnan velocities. Since M r 2 ~^~ d 
is a contact term oc 5(r) for k — 0, one has d t (R 2 ^ d ) oc V(r;0;t). The latter is zero 
in the smooth case so that (R~ d ) is a true integral of motion. In the nonsmooth case, 
(R 2 ~^ d ) oc t 1 ~ d /( 2 ~5) anc [ j s no t conserved due to a nonzero probability density to find two 
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particles at the same place even when they started apart. 

As stated, the result (55) holds when the molecular diffusivity is turned off in the velocity 
ensemble with no viscous cutoff, i.e. for vanishing Schmidt number Sc = vj k, where v is the 
viscosity defined in the Kraichnan model as D\°cf-. The same result holds also when v and 
k are turned off at the same time with Sc = O(l), provided the initial distance r is taken 
to zero only afterwards (E and Vanden Eijnden, 2000b). This confirms that the explosive 
separation of close trajectories persists for finite Reynolds numbers as long as their initial 
distance is not too small, as anticipated by Bernard et al. (1998). 

D. Two-particle dispersion in a compressible flow 

Discussing the particle dispersion in incompressible fluids and exposing the different 
mechanisms of particle separation, we paid little attention to the detailed geometry of the 
flows, severely restricted by the incompressibility. The presence of compressibility allows 
for more flexible flow geometries with regions of ongoing compression effectively trapping 
particles for long times and counteracting their tendency to separate. To expose this ef- 
fect and gauge its relative importance for smooth and nonsmooth flows, we start from the 
simplest case of a time-independent Id flow x = v(x) . In Id, any velocity is potential: 
v(x) = —d x (f)(x) , and the flow is the steepest descent in the landscape defined by the po- 
tential 0. The particles are trapped in the intervals where the velocity has a constant sign 
and they converge to the fixed points with lower value of at the ends of those intervals. In 
the regions where d x v is negative, nearby trajectories are compressed together. If the flow 
is smooth the trajectories take an infinite time to arrive at the fixed points (the particles 
might also escape to infinity in a finite time). Let us consider now a nonsmooth version of 
the velocity, e.g. a Brownian path with Holder exponent 1/2. At variance with the smooth 
case, the solutions will take a finite time to reach the fixed points at the ends of the trap- 
ping intervals and will stick to them at subsequent times, as in the example of the equation 
X = \x — Xq | 1/2 . The nonsmoothness of the velocity clearly amplifies the trapping effects 
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leading to the convergence of the trajectories. A time-dependence of the velocity changes 
somewhat the picture. The trapping regions, as defined for the static case, start wandering 
and they do not enslave the solutions which may cross their boundaries. Still, the regions 
of ongoing compression effectively trap the fluid particles for long time intervals. Whether 
the tendency of the particles to separate or the trapping effects win is a matter of detailed 
characteristics of the flow. 

In higher dimensions, the behavior of potential flows is very similar to the Id case, with 
trapping totally dominating in the time-independent case, its effects being magnified by 
the nonsmoothness of the velocity and blurred by the time- dependence. The traps might 
of course have a more complicated geometry. Moreover, we might have both solenoidal 
and potential components in the velocity. The dominant tendency for the incompressible 
component is to separate the trajectories, as we discussed in the previous Sections. On the 
other hand, the potential component enhances trapping in the compressed regions. The 
net result of the interplay between the two components depends on their relative strength, 
spatial smoothness and temporal rate of change. 

Let us consider first a smooth compressible flow with a homogeneous and stationary er- 
godic statistics. Similarly to the incompressible case discussed in Sect. II.B.l, the stretching- 
contraction variables pi, i = 1, . . .,d, behave asymptotically as t\i with the PDF of large 

deviations Xi = pi/t — \ determined by an entropy function H(xi, . . . , xj). The asymptotic 

d 

growth rate of the fluid volume is given by the sum of the Lyapunov exponents s — \- 

i=i 

Note that density fluctuations do not grow in a statistically steady compressible flow because 
the pressure provides feedback from the density to the velocity field. That means that s 
vanishes even though the pi variables fluctuate. However, to model the growth of density 
fluctuations in the intermediate regime, one can consider an idealized model with a steady 
velocity statistics having nonzero s. This quantity has the interpretation of the opposite 
of the entropy production rate, see Section III. A. 4 below, and it is necessarily < (Ruelle, 
1997). Let us give here the argument due to Balkovsky et al. (1999a) which goes as follows. 
In any statistically homogeneous flow, incompressible or compressible, the distribution of 
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particle displacements is independent of their initial position and so is the distribution of 
the evolution matrix Wij(t;r) = d R 1 (t; r) / dr j . Since the total volume V (assumed finite 
in this argument) is conserved, the average (det W) is equal to unity for all times and ini- 
tial positions although the determinant fluctuates in the compressible case. The average of 
det W = e S a is dominated at long times by the saddle-point x* giving the maximum of 
J2(K + Xi) — H (x) , which has to vanish to conform with the total volume conservation. Since 
J2%i — H{x) is concave and vanishes at x — 0, its maximum value has to be non-negative. 
We conclude that the sum of the Lyapunov exponents is non-positive. The physics behind 
this result is transparent: there are more Lagrangian particles in the contracting regions, 
leading to negative average gradients in the Lagrangian frame. Indeed, the volume growth 
rate tends at large times to the Lagrangian average of the trace of the strain a = Vv: 

- t (\ndet (W(r;t))) — J tr a(R(t; r), t) % = J tva(R,t) - J*^ . (56) 

The Lagrangian average generally coincides with the Eulerian one / dr tr a(t, r)/V, only in 
the incompressible case (where it is zero). For compressible flow, the integrals in (56) vanish 
at the initial time (when we set the initial conditions for the Lagrangian trajectories so that 
the measure was uniform back then) . The regions of ongoing compression with negative tr a 
acquire higher weight in the average in (56) than the expanding ones. Negative values of 
tr a suppress stretching and enhance trapping and that is the simple reason for the volume 
growth rate to be generally negative. Note that, were the trajectory R(r; t) defined by its 
final (rather than initial) position, the sign of the average strain trace would be positive. Let 
us stress again the essential difference between the Eulerian and the Lagrangian averages 
in the compressible case: an Eulerian average is uniform over space, while in a Lagrangian 
average every trajectory comes with its own weight determined by the local rate of volume 
change. For the corresponding effects on the single-particle transport, the interested reader 
is referred to (Vergassola and Avellaneda, 1997). 

In the particular case of a short-correlated strain one can take t in (56) larger than the cor- 
relation time r of the strain, yet small enough to allow for the expansion det (W(r; t)) _1 ~ 
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1 — f tra(r,t') dt' so that (56) becomes equal to — J*(tr a(r, t) tra(r, t')) dt' . More formally, 

o 

let us introduce the compressible generalization of the Kraichnan ensemble for smooth ve- 
locities. Their (non-constant part of the) pair correlation function is defined as 

d ij (r) = D x [(d + 1 - 2p) 5 ij r 2 + 2(pd - 1) rV] + o(r 2 ) , (57) 

compare to (29). The degree of compressibility p = ((Vji> l ) 2 )/((Vji> J ) 2 ) is between and 1 
for the isotropic case at hand, with the the two extrema corresponding to the incompressible 
and the potential cases. The corresponding strain matrix a = Vv has the Eulerian mean 
equal to zero and 2-point function 

( <Jij(t) a u {t') ) = 25{t- t') D 1 [(d+1- 2p) 5 lk 5 j£ + (pd - 1)(%<W + S u S jk )] . (58) 

t 

The volume growth rate — / (au(t) Gjj{t')) dt 1 is thus strictly negative, in agreement with 

o 

the general discussion, and equal to —pDid(d — l)(d + 2) if we set / 5(t)dt = 1/2. 
The same result is obtained more systematically by considering the Ito stochastic equation 
dW = adtW for the evolution matrix and applying the Ito formula to lndet(J¥), see 
Appendix. One may identify the generator of the process W(t) and proceed as for the 
incompressible case calculating the PDF V(pi, ■ ■ ■ , Pd',t). It takes again the large deviation 
form (21), with the entropy function and the Lyapunov exponents given by 



H ^ ~ 4Di(*f P (d-2)) [S x i + P (d-i)(d+2) ( It, Xi ) } ' ( 59 ) 

i=l i=l 

Xi = D x [d(d - 2i + 1) - 2p (d + (d - 2)i)] , (60) 



to be compared to (32). Note how the form (59) of the entropy imposes the condition 
= in the incompressible limit. The interparticle distance R(t) has the lognormal 
distribution (23) with X = X 1 = D^d- l)(d- Ap) and A = 2D 1 (d-l)(l + 2p). Explicitly, 
t- 1 ln(i? n ) oc n[n + d + 2p{n - 2)\ (Chertkov et ai, 1998). The quantity r(4p-W(Wp) [ s 
thus statistically conserved. The highest Lyapunov exponent A becomes negative when the 
degree of compressibility is larger than d/A (Le Jan, 1985; Chertkov et ai, 1998). Low-order 



40 



moments of R, including its logarithm, would then decrease while high-order moments would 
grow with time. 

It is instructive to decompose the strain into its "incompressible and compressible parts" 
(Tij — ^Sijtrcr and ^5ijtra. From the equality A« = (an), see (21), it follows that the 
Lyapunov exponents of the incompressible part (having A > 0) get uniformly shifted down by 
the Lagrangian average of tra/d. In Id, where the compressibility is maximal 4 , A < 0. The 
lowering of the Lyapunov exponents when p grows clearly signals the increase of trapping. 
The regime with p > d/4 , with all the Lyapunov exponents becoming negative, is the one 
where trapping effects dominate. The dramatic consequences for the scalar fields advected 
by such flow will be discussed in Sect. III.B.l. 

As it was clear from the Id example, we should expect even stronger effects of com- 
pressibility in nonsmooth velocity fields, with an increased tendency for the fluid particles 
to aggregate in finite time. This has, indeed, been shown to occur when the velocity has a 
short correlation time, i.e. for the nonsmooth version of the compressible Kraichnan ensem- 
ble (Gaw§dzki and Vergassola, 2000). The expression (46) for the Fourier transform of the 
2-point correlation function is easily modified. The functional dependence on k 2 remains 
the same, the solenoidal projector is simply multiplied by 1 — p and one adds to it the 
compressible longitudinal component p(d — l)k l ki jk 2 . This gives for the non-constant part 
of the 2-point function 

tfJ(r) = D ± [(d - 1 + £ - p£) 5 ij + £(pd - 1) rW~ 2 ] . (61) 

For £ = 2, (61) reproduces (57) without the o(r 2 ) corrections. Most of the results discussed 
in Sect. II. C for the incompressible version of the model still go through, including the con- 
struction of the Markov process describing the noisy trajectories and the heat-kernel form 
of the joint PDF of two particle positions. The restriction M of M 2 to the homogeneous 
and isotropic sector, whose heat kernel gives the PDF V(r; R; t) of the distance between 

4 One-dimensional results are recovered from our formulae by taking p = 1 and D\ oc l/(d — 1). 
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two particles, takes now the form 



M 



j,a— 1 — 7 yd—Y 



d r , (62) 



where 7 = £>£(d + £)/(l + p£). As found in (Gaw§dzki and Vergassola, 2000), see also 
(Le Jan and Raimond, 1999; E and Vanden Eijnden, 2000b), depending on the smoothness 
exponent £ and the degree of compressibility p, two different regimes arise in the limit k — > 0. 

For weak compressibility p < p c = |j, the situation is very much the same as for the 
incompressible case and 

lim V(r;R;t) 

c*- I w/o c\ ex P ~~ const. — — — . (63) 



|(|(d-7)/(2-0 

We still have an explosive separation of trajectories but, in comparison to the incompressible 
situation, the power prefactor Rr~< with 7 oc p suppresses large separations and enhances 
small ones. When p crosses p c = d+ 2 ^~ 2 , the particle-touching event R = becomes 
recurrent for the Markov process describing the distance between the two particles (Le Jan 
and Raimond, 1999). In other words, for p c < p < p c a pair of Lagrangian trajectories 
returns infinitely often to a near touch, a clear sign of increased trapping. 

When p crosses p c , the singularity at R — of the right hand side of (63) becomes 
non-integrable and a different limit is realized. Indeed, in this regime, 

lim V(r;R;t) = V re9 (r; R;t) + p(r;t)5(R), (64) 

with the regular part V reg tending to zero and p approaching unity when r — > 0. This re- 
produces for i] = the result (54), always holding when the viscous cutoff r\ > smoothes the 
velocity realizations. In other words, even though the velocity is nonsmooth, the Lagrangian 
trajectories in a fixed velocity field are determined by their initial positions. Moreover, the 
contact term in (64) signals that trajectories starting at a finite distance r collapse to zero 
distance and stay together with a positive probability growing with time (to unity if the 
integral scale L — 00). The strongly compressible regime p > p c is clearly dominated by 
trapping effects leading to the aggregation of fluid particles, see Fig. 2. The same results 
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hold if we turn off the diffusivity k and the viscosity v at the same time, with the notable 
exception of the intermediate regime p c < p < p c . In this interval, if the Schmidt num- 
ber Sc diverges fast enough during the limiting process, the resulting PDF of the distance 
takes the form (64) rather than that of (55) arising for bounded Sc (E and Vanden Eijnden, 
2000b). Sufficiently high Schmidt numbers thus lead to the particle aggregation in this case. 
Note that in the limit of smooth velocities £ — > 2, the intermediate interval shrinks to the 
point p — d/4 where the highest Lyapunov exponent crosses zero. 

As was mentioned, the aggregation of fluid particles can take place only as a transient 
process. The back reaction of the density on the flow eventually stops the growth of the 
density fluctuations. The transient trapping should, however, play a role in the creation of 
the shocklet structures observed in high Mach number compressible flows (Zeman, 1990). 
There is another important physical situation that may be modeled by a smooth compressible 
random flow with a nonzero sum of the Lyapunov exponents. Let us consider a small inertial 
particle of density p and radius a in a fluid of density p . Its movement may be approximated 
by that of a Lagrangian particle in an effective velocity field provided that a 2 jv is much 
smaller than the velocity time scale in the Lagrangian frame. The inertial difference between 
the effective velocity v of the particle and the fluid velocity u(r,t) is proportional to the 
local acceleration: v = u + {(5 — l)r s du/dt, where (3 = 3p/(p + 2p ) and r s = a 2 /3z//9 is 
the Stokes time. Considering such particles distributed in the volume, one may define the 
velocity field v(r, t) , whose divergence oc V[(tt • V)w] does not vanish even if the fluid flow 
is incompressible. As discussed above, this leads to a negative volume growth rate and the 
clustering of the particles (Balkovsky et ai, 2001). 

E. Multiparticle dynamics, statistical conservation laws and breakdown of scale invariance 

This subsection is a highlight of the review. We describe here the time-dependent statis- 
tics of multiparticle configurations with the emphasis on conservation laws of turbulent trans- 
port. As we have seen in the previous subsections, the two-particle statistics is characterized 
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by a simple behavior of the single separation vector. In nonsmooth velocities, the length of 
the vector grows by a power law, while the initial separation is forgotten and there are no 
statistical integrals of motion. In contrast, the many-particle evolution exhibits non-trivial 
statistical conservation laws that involve geometry and are proportional to positive powers 
of the distances. The distance growth is balanced by the decrease of the shape fluctuations 
in those integrals. The existence of multiparticle conservation laws indicates the presence of 
a long-time memory and is a reflection of the coupling among the particles due to the simple 
fact that they all are in the same velocity field. The conserved quantities may be easily built 
for the limiting cases. For very irregular velocities, the fluid particles undergo independent 
Brownian motions and the interparticle distances grow as (R^ m (t)) = R^ m (0)+Dt. Here, ex- 
amples of statistical integrals of motion are (R 2 nm —R 2 pr ) and (2(d+2)R 2 im R^ r —d(R^ im +Rp r )) : 
and an infinity of similarly built harmonic polynomials where all the powers of t cancel out. 
Another example is the infinite-dimensional case, where the interparticle distances do not 
fluctuate. The two-particle law (38), R]^(t) — Rn~^*(0) oc t, implies then that the ex- 
pectation of any function of R]^ — R 1 ^ does not change with time. A final example is 
provided by smooth velocities, where the particle separations at long times become aligned 
with the eigendirections of the largest Lyapunov exponent of the evolution matrix W(t) 
defined in (16). All the interparticle distances R nm will then grow exponentially and their 
ratios R nm /Rki do not change. Away from the degenerate limiting cases, the conserved 
quantities continue to exist, yet they cannot be constructed so easily and they depend on 
the number of particles and their configuration geometry. The very existence of conserved 
quantities is natural; what is generally nontrivial is their precise form and their scaling. 
The intricate statistical conservation laws of multiparticle dynamics were first discovered for 
the Kraichnan velocities. That came as a surprise since the Kraichnan velocity ensemble is 
Gaussian and time-decorrelated, with no structure built in, except for the spatial scaling in 
the inertial range. The discovery has led to a new qualitative and quantitative understand- 
ing of intermittency, as we shall discuss in detail in Sect. III.C.l. Even more importantly, it 
has pointed to aspects of the multiparticle evolution that seem both present and relevant in 
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generic turbulent flows. Note that those aspects are missed by simple stochastic processes 
commonly used in numerical Lagrangian models. There is, for example, a long tradition 
to take for each trajectory the time integral of a (/-dimensional Brownian motion (whose 
variance is oc t 3 as in the Richardson law) or an Ornstein-Uhlenbeck process. Such models, 
however, cannot capture correctly the subtle features of the iV-particle dynamics such as 
the statistical conservation laws. 

1. Absolute and relative evolution of particles 

As for many-body problems in other branches of physics, e.g. in kinetic theory or in 
quantum mechanics, multiparticle dynamics brings about new aspects due to cooperative 
effects. In turbulence, such effects are mediated by the velocity fluctuations with long space- 
correlations. Consider the joint PDF's of the equal-time positions R = (R 1 , . . . , R N ) of A" 
fluid trajectories 

N 

(Y[p(r n ,s; R n ,t\v)) = V N (r; R; t — s), (65) 

71=1 

with the average over the velocity ensemble, see Fig. 3. More generally, one may study the 
different-time versions of (65). Such PDF's, called multiparticle Green functions, account 
for the overall statistics of the many-particle systems. 

For statistically homogeneous velocities, it is convenient to separate the absolute motion 
of the particles from the relative one, as in other many-body problems with spatial homo- 
geneity. For a single particle, there is nothing but the absolute motion which is diffusive 
at times longer that the Lagrangian correlation time (Sect. II. A). For A" particles, we may 
define the absolute motion as that of the mean position R = J^Rn/N. When the parti- 

n 

cles separate beyond the velocity correlation length, they are essentially independent. The 
absolute motion is then diffusive with the diffusivity N times smaller than that of a single 
particle. The relative motion of A" particles may be described by the versions of the joint 
PDF's (65) averaged over rigid translations: 
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V N {r;E;t) = JV N (r;E+p;,t)dp, (66) 

where p = (p, . . . , p). The PDF in (66) describes the distribution of the particle separations 
Rnm = Rn — Rm or of the relative positions R = (R 1 — R, . . . , .R N — R). 

As for two particles, we expect that when k — > the multiparticle Green functions 
V N tend to (possibly distributional) limits that we shall denote by the same symbol. The 
limiting PDF's are again expected to show a different short-distance behavior for smooth 
and nonsmooth velocities. For smooth velocities, the existence of deterministic trajectories 
leads for k = to the collapse property 

lim V N (r; R, t) = V N _^ R, t) S(R N ^ - R N ), (67) 

N N-1 

where R = (R 1 , . . . , R N _ 1 ) and similarly for the relative PDF's. If all the distances among 
the particles are much smaller than the viscous cutoff, the velocity differences are approxi- 
mated by linear expressions and 

r N 

V N (r;R;t) = (l[6(R n + p- W(t) r n ) ) dp . (68) 

J n=l 

The evolution matrix W(t) was defined in (16) and the above PDF's clearly depend only on 
its statistics which has been discussed in Sect. II. B. 

2. Multiparticle motion in Kraichnan velocities 

The great advantage of the Kraichnan model is that the statistical Lagrangian integrals 
of motion can be found as zero modes of explicit evolution operators. Indeed, the crucial 
simplification lies in the Markov character of the Lagrangian trajectories due to the velocity 
time decorrelation. In other words, the processes R(t) and R(t) are Markovian and the 
multiparticle Green functions V N and V N give, for fixed N, their transition probabilities. The 
process R(t) is characterized by its second-order differential generator M N , whose explicit 
form may be deduced by a straightforward generalization of the path integral representation 
(52) to N particles. The PDF V N (r; R, t) = el'- s l A/1 iv(r; R) with 
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N N 

M N = E Dij (rnm)V ri V rL + V '„- (69) 

n,m=l n=l 

For the relative process R(t), the operator TW^ should be replaced by its translation- 
invariant version 

Ms = ~ E (>>nj + 2«(5«) V r , V^, (70) 

n<m 

with cf J related to D' lj by (47). Note the multibody structure of the operators in (69) and 
(70). The limiting PDF's obtained for k — > define the heat kernels of the k — version of 
the operators that are singular elliptic and require some care in handling (Hakulinen, 2000). 

As we have seen previously, the Kraichnan ensemble may be used to model both smooth 
and Holder continuous velocities. In the first case, one keeps the viscous cutoff rj in the 
two-point correlation (46) with the result that d^{r) = 0(r 2 ) for r <C r\ as in (29), or one 
sets £ = 2 in (48). The latter is equivalent to the approximation (68) with W(t) becoming 
a diffusion process on the group SL(d) of unimodular matrices, with an explicitly known 
generator, as discussed in Sect. II.B.2(i). The right hand side of (68) may then be studied 
by using the representation theory (Shraiman and Siggia, 1995 and 1996; Bernard et ai, 
1998), see also Sect. II. E. 5 below. It exhibits the collapse property (67). 

From the form (70) of the generator of the process R(t) in the Kraichnan model, we infer 
that N fluid particles undergo an effective diffusion with the diffusivity depending on the 
interparticle distances. In the inertial interval and for a small molecular diffusivity k, the 
effective diffusivity scales as the power £ of the interparticle distances. Comparing to the 
standard diffusion with constant diffusivity, it is intuitively clear that the particles spend 
longer time together when they are close and separate faster when they become distant. Both 
tendencies may coexist and dominate the motion of different clusters of particles. It remains 
to find a more analytic and quantitative way to capture those behaviors. The effective short- 
distance attraction that slows down the separation of close particles is a robust collective 
phenomenon expected to be present also in time-correlated and non-Gaussian velocity fields. 
We believe that it is responsible for the intermittency of scalar fields transported by high 
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Reynolds number flows, as it will be discussed in the second part of the review. 

As for a single particle, the absolute motion of N particles is dominated by velocity 
fluctuations on scales of order L. In contrast, the relative motion within the inertial range 
is approximately independent of the velocity cutoffs and it is convenient to take directly the 
scaling limit rj = and L = oo. We shall also set the molecular diffusivity to zero. In these 
limits, M. N has the dimension length^ -2 , implying that time scales as length 2 ~^ and 

V N (Xr; R- t) = A _(JV_1)d V N (r; X^R; X^H) . (71) 

The relative motion of iV fluid particles may be tested by tracing the time evolution of the 
Lagrangian averages 

( / (R(t)) ) = / f(R) V N (r; R, t) dR (72) 

of translation-invariant functions / of the simultaneous particle positions. Think about the 
evolution of iV fluid particles as that of a discrete cloud of marked points in physical space. 
There are two elements in the evolution of the cloud: the growth of its size and the change of 
its shape. We shall define the overall size of the cloud as R = (^ Rnm) 1 ^ 2 an d its "shape" 
as R = R/ R. For example, three particles form a triangle in space (with labeled vertices) 
and the notion of shape that we are using includes the rotational degrees of freedom of the 
triangle. The growth of the size of the cloud might be studied by looking at the Lagrangian 
average of the positive powers -R c . More generally, let / be a scaling function of dimension 
C, i.e. such that f(\R) = \** f(R)- The change of variables -R h- > t^R, the relation (71) 
and the scaling property of / allow to trade the Lagrangian PDF V N in (72) for that at 
unit time V N (t~~r; R, 1). As for two points, the limit of V N when the initial points 
approach each other is nonsingular for nonsmooth velocities and we infer that 

(f(R(t)) ) = t& J f(R) V N (0; R, 1) dR' + o(t&). (73) 

In particular, we obtain the iV-particle generalization of the Richardson- type behavior (39): 

( R(t)^ ) oc t~ . Hence, in the Kraichnan model the size of the cloud of Lagrangian points 

i 

grows superdiffusively cx t 2 ^. What about its shape? 
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3. Zero modes and slow modes 



In order to test the evolution of the shape of the cloud one might compare the Lagrangian 
averages of different scaling functions. The relation (73) suggests that at large times they all 
scale dimensionally as t~ . Actually, all do but those for which the integral in (73) vanishes. 
The latter scaling functions, whose evolution violates the dimensional prediction, may thus 
be better suited for testing the evolution of the shape of the cloud. Do such functions 
exist? Suppose that / is a scaling function of non-negative dimension ( annihilated by Ai N , 
i.e. such that M. N f = 0. Its Lagrangian average, rather than obeying the dimensional law, 
does not change in time: (R(t)) ^ = f(r). Indeed, d t V N (r; R; t) = M N V N (r; R; t). 
Therefore, the time-derivative of the right hand side of (72) vanishes since it brings down the 
(Hermitian) operator A4 N acting on /. Thus, the zero modes of M. N are conserved in mean 
by the Lagrangian evolution. The importance of such conserved modes for the transport 
properties by 5-correlated velocities has been recognized independently by Shraiman and 
Siggia (1995), Chertkov et at. (1995b), and Gaw§dzki and Kupiainen (1995 and 1996). 

The above mechanism may be easily generalized (Bernard et al, 1998). Suppose that 
fk is a zero mode of the (k + l) st power of M. N (but not of a lower one), with scaling 
dimension ( + (2 — Then, its Lagrangian average is a polynomial of degree k in time 
since its (k+l) st time derivative vanishes. Its temporal growth is slower than the dimensional 
prediction t^ +k if ( > so that the integral coefficient in (73) must vanish. We shall call 
such scaling functions slow modes. The slow modes may be organized into "towers" with 
the zero modes at the bottom 5 . One descends down the tower by applying the operator M N 
which lowers the scaling dimension by (2 — £). The zero and the slow modes are natural 
candidates for probes of the shape evolution of the Lagrangian cloud. There is an important 
general feature of those modes due to the multibody structure of the operators: the zero 
modes of A1 JV _ 1 are also zero modes of A4 N and the same for the slow modes. Only those 

5 Note that (M N ) k fk is a zero mode of scaling dimension £. 
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modes that depend non-trivially on all the positions of the N points may of course give 
new information on the iV-particle evolution which cannot be read from the evolution of a 
smaller number of particles. We shall call such zero and slow modes irreducible. 

To get convinced that zero and slow modes do exist, let us first consider the limiting 
case £ — > of very rough velocity fields. In this limit d tj (r) reduces to Di(d — l)5' lj , see 
(48), and the operator Ai N becomes proportional to V 2 , the (iV<i)-dimensional Laplacian 
restricted to the translation-invariant sector. The relative motion of the particles becomes 
a pure diffusion. If R denotes the size-of-the-cloud variable then 

V 2 = R- d ^ +1 d R R d ^- 1 d R + R- 2 V 2 , (74) 

- — -2 

where d N = (N—l)d and V is the angular Laplacian on the (d N — l)-dimensional unit sphere 
of shapes R. The spectrum of the latter may be analyzed using the properties of the rotation 
group. Its eigenvalues are + d N — 2), where j — 0, 1, . . . is the angular momentum. 

The functions / J)0 = Ri<f>j(R), where <f>j is an angular momentum j eigenfunction, are zero 
modes of the Laplacian with scaling dimension j. The contributions coming from the radial 
and the angular parts in (74) indeed cancel out. The polynomials fj jk = R 2k fj :0 form the 
corresponding (infinite) tower of slow modes. All the scaling zero and slow modes of the 
Laplacian are of that form. The mechanism behind the special behavior of their Lagrangian 
averages is that mentioned at the beginning of the section. Beside the constant, the simplest 
zero mode has the form of the difference R 2 2 — -R 2 3 . Both terms are slow modes in the tower 
of the constant zero mode and their Lagrangian averages grow linearly in time with the same 
leading coefficient. Their difference is thus constant. A similar mechanism stands behind 
the next example, the difference 2(d + 2)R 2 2 R 2 4 — d(Rf 2 + -R34), whose Lagrangian average 
is conserved due to the cancellation of linear and quadratic terms in time, and so on. 

It was argued by Bernard et al. (1998) that the slow modes exist also for general £ and 
show up in the asymptotic behavior of the multiparticle PDF's when the initial points get 
close: 
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oo 

V N (Xr; E;t) = EE f ak{ r) g a , k (R,t), (75) 

a k=0 

for small A. The first sum is over the zero modes f a;0 = f a with scaling dimensions ( a , 
while higher k give the corresponding towers of slow modes. The functions in (75) may be 
normalized so that / a , fc _i = M N f a , k and g a , k+1 = M N g a ,k = d t g a ,k- The leading term in the 
expansion comes from the constant zero mode /o,o = 1- The corresponding g ,o coincides 
with the PDF of N initially overlapping particles. The asymptotic expansion (75) is easy to 
establish for vanishing £ and for N = 2 with arbitrary £ . In the general case, it has been 
obtained under some plausible, but yet unproven, regularity assumptions. Note that, due 
to (71), the expansion (75) describes also the asymptotics of the multiparticle PDF's when 
the final points get far apart and the times become large. The use of (75) allows to extract 
the complete long-time asymptotics of the Lagrangian averages: 

oo „ 

( / (R(t)) ) = E E U{r) / f(E) 9a*(R, 1) d# , (76) 

which is a detailed refinement of (73), corresponding just to the first term. Note that the 
pure polynomial-in-time behavior of the Lagrangian averages of the slow modes implies 
partial orthogonality relations between the slow modes and the g modes. 

4. Shape evolution 

The qualitative mechanism behind the preservation of the Lagrangian average of the zero 
modes is the compensation between its increase due to the size growth and its depletion due 
to the shape evolution. The size and the shape dynamics are mixed in the expansion (75) 
that, together with (71), describes the long-time long-distance relative evolution of the 
Lagrangian cloud. To get more insight into the cooperative behavior of the particles and 
the geometry of their configurations, it is useful to separate the shape evolution following 
Gat and Zeitak (1998), see also (Arad and Procaccia, 2000). The general idea is to trade 
time in the relative iV-particle evolution R(t) for the size variable R(t). This may be done 
as follows. Let us start with iV particles in a configuration of size r and shape r. Denote 
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by R(R) the shape of the particle configuration the first time its size reaches R > r. 
Varying R, one obtains a description of the evolution of the shape of the particle cloud 
with its size (which may be discontinuous if the size does not grow at all moments of 
time). For scale invariant velocity fields, the PDF of the shapes R(R), i.e. of the first 
passages of R(t) through the sphere of size R, depends only on the ratio R/r. We shall 
denote it by P N (r;R; r/R) . The shape evolution R(R) is still a Markov process: in order 
to compute the probability of the first passage through the sphere of size R, one may 
condition with respect to the first passage through a sphere of an intermediate size. As a 
result, the PDF's P N obey a semigroup Chapman-Kolmogorov relation. As observed by 
Gat and Zeitak (1998), the eigenmode expansion of the (generally non-Hermitian) Markov 
semigroup P N (X) involves the zero modes f a : 

PAtRA) = E X(a fad) K(R). (77) 

a 

The formal reason is as follows. The statistics of the first passage through a given surface may 
be obtained by imposing the Dirichlet boundary condition on the surface in the differential 
generator of the process. For the case at hand, if G N (r , R) denotes the kernel of the inverse 
of —M AT and G Dir refers to its version with the Dirichlet conditions at R—l, then 

N 

P N (t,R;X) = \ E [v/][^(i4 m )V C(Ar,E)], (78) 

n,m=l 

with the derivatives taken on the sphere of unit radius. The potential theory relation (78) 
expresses the simple fact that the probability of a first passage through a given surface is 
the normal component of the probability current (the expression in second parenthesis on 
the right hand side). On the other hand, by integrating the asymptotic expansion (75), one 
obtains the expansion 

oo 

G N (Xr,R) = fv N (Xr;R,t)dt = J2^ a fa(r)g a (R), (79) 
o 

oo 

where g a = J g a ^ dt. Note that the slow modes do not appear since g a ,k+i is the time 
o 

derivative of g a ^, that vanishes at the boundary of the integration interval for nonzero 
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R. The Dirichlet boundary condition at it! = 1 should not affect the dependence on the 
asymptotically close positions Xr. An expansion analogous to (79) should then hold for 
G° ir with the same zero modes and modified functions g. The potential theory formula 
(78) together with the Dirichlet version of (79) give (77). 

Let us now consider the "shape only" version of the Lagrangian average (72). Substitut- 
ing the expansion (77), we obtain 

(/ (R(R))) = E ( R M~ Ca fair.) I /(I) K(R) dR. (80) 

a J 

The interpretation of (80) is simple: when the size increases, the average of a generic function 
of the shape relaxes as a combination of negative powers of R to a constant, with the zero 
modes f a (z) giving the modes of relaxation (Gat and Zeitak, 1998). On the other hand, 
due to the orthogonality of the left and right eigenfunctions of the semigroup P/v(A), the 

shape averages of the zero modes decay: (^f a (£L(R))^ = (R/r)^ a f a (r_). Since the size 

i 

grows as t 2 -«, this quantitatively illustrates the decrease of the shape average responsible 
for the conservation in time of the Lagrangian average of the zero modes. A vivid and 
explicit example of the compensation is provided by the case of three particles. The contour 
lines of the relevant zero mode as a function of the shape of the triangle are shown in Fig. 4. 
The function tends to decrease for configurations where all the interparticle distances are 
comparable. It is then clear that the decrease in the shape average is simply due to particle 
evolving toward symmetrical configurations with aspect ratios of order unity (Pumir, 1998; 
Celani and Vergassola, 2001). 

The relative motion of the particles in the limit £ — > becomes a pure diffusion. It 
is then easy to see that the zero modes of the Laplacian play indeed the role of relaxation 
modes of the £ = shape evolution. Pure diffusion is the classical case of potential theory: 
G N (r,R) oc \r — R\~ d N + 2 is the potential induced at r by a unit charge placed at R 
(the absolute value is taken in the sense of the size variable). The Dirichlet version G^ ir 
corresponds to the potential of a unit charge inside a grounded conducting sphere and it 
is obtained by the image charge method: G^ ir (r, i?) = G N (r,R) — R~ d N + 2 G N (r,R/R 2 ). 
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The potential theory formula (78) gives then for the shape-to-shape transition probability 
-Pjv(£; R'j A) oc (1 — A 2 ) \R— \r_\~ dN , with the proportionality constant equal to the inverse 
volume of the sphere. On the other hand, calculating G^ ir in each angular momentum 
sector by (74), one easily describes the generator of the Markov semigroup P N {\). It is 
the pseudo-differential Hermitian operator with the same eigenf unctions 4>j as the angular 
Laplacian V but with eigenvalues —j. The expansion (77) with h a = f a follows (recall 
that the functions RP<f>j form the zero modes of the Laplacian). The Markov process R(R) 
lives on distributional realizations (and not on continuous ones). 

It is instructive to compare the shape dynamics of the Lagrangian cloud to the imaginary- 
time evolution of a quantum mechanical many-particle system governed by the Hamiltonian 

2 

H N = J2 + V(r nm ). The (Hermitian) imaginary-time evolution operators e" 

n n<m 

decompose in the translation-invariant sector as J2 e ~ tEN,a |VVa)(VvJ • The ground state 

a 

energy is E N and the sum is replaced by an integral for the continuous part of the spectrum. 
An attractive potential between the particles may lead to the creation of bound states at 
the bottom of the spectrum of H N . Breaking the system into subsystems of Ni particles 
by removing the potential coupling between them would then raise the ground state energy: 
E N o < S-^jv o- A ver y similar phenomenon occurs in the stochastic shape dynamics. Con- 
sider indeed an even number of particles and denote by ( N0 the lowest scaling dimension 
of the irreducible zero mode invariant under translations, rotations and reflections. For two 
particles there is no invariant irreducible zero mode and its role is played by the first slow 
mode oc r 2 ~^ and ( 20 = 2 — £. Suppose now that we break the system into subsystems with 
an even number iVj of particles by removing in A4 N the appropriate terms rf(r nm )V rn V rm 
coupling the subsystems, see (70). For Ni > 4, the irreducible zero mode for the broken 
system factorizes into the product of such modes for the subsystems. If iVj = 2 for some 
i, the factorization still holds modulo terms dependent on less variables. In any case, the 
scaling dimensions simply add up. The crucial observation, confirmed by perturbative and 
numerical analyses discussed below, is that the minimal dimension of the irreducible zero 



54 



modes is raised: ( N < Yl C N - o ■ m particular, is smaller than ^(2 — £). One even 
expects that ( N0 i s a concave function of (even) N and that for N ^> d its dependence on 
iV saturates, see Sects. III. C. 2, III. D. 2, III. F. By analogy with the multibody quantum me- 
chanics, we may say that the irreducible zero modes are bound states of the shape evolution 
of the Lagrangian cloud. The effect is at the root of the intermittency of a passive scalar 
advected by nonsmooth Kraichnan velocities, as we shall see in Sect. III. C.l. It is a coop- 
erative phenomenon exhibiting a short-distance attraction of close Lagrangian trajectories 
superposed on the overall repulsion of the trajectories. There are indications that similar 
bound states of the shape evolution persist in more realistic flows and that they are still 
responsible for the scalar intermittency, see Sect. III.D.2 and (Celani and Vergassola, 2001). 

5. Perturbative schemes 

The incompressible Kraichnan model has two parameters d G [2, oo) and £ e [0, 2]. It is 
then natural to ask if the problem is simplified at their limiting values and if perturbative 
methods might be used to get the zero modes near the limits. No significant simplification 
has been recognized for d — 2 at arbitrary £. The other limits do allow for a perturbative 
treatment since the particle interaction is weak and the anomalous scaling disappears there. 
The perturbation theory is essentially regular for £ ^ and d — > oo. Conversely, the 
perturbation theory for £ — > 2 is singular for two reasons. First, the advection by a smooth 
velocity field preserves the configurations with the particles aligned on a straight line. A 
small roughness of the velocity has little effect on the particle motion almost everywhere 
but for quasi-collinear geometries. A separate treatment of those regions and a matching 
with the regular perturbation expansion for a general geometry is thus needed. Second, for 
almost smooth velocities, close particles separate very slowly and their collective behavior 
is masked by this effect which leads to an accumulation of zero modes with very close 
scaling dimensions. We shall start by the more regular cases of small £ and large d. The 
scaling of the irreducible four-point zero mode with the lowest dimension was first calculated 
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to the linear order in £ by Gaw§dzki and Kupiainen (1995) by a version of degenerate 
Rayleigh-Schrodinger perturbation theory. In parallel, a similar calculation in the linear 
order in 1/d was performed by Chertkov et al. (1995b). Bernard et al. (1996) streamlined 
the small £ analysis and generalized it to any even order, following a similar generalization 
by Chertkov and Falkovich (1996) for the 1/d expansion. We sketch here the main lines of 
those calculations. 

As we discussed in Section II. E. 3, the operator M N is reduced to the Laplacian (74) for 
£ = 0. The zero modes of the Laplacian depend on the size of the particle configuration as 
E? and on its shape as the eigenfunctions of V with angular momentum j. The zero modes 
invariant under d- dimensional translations, rotations and reflections can be reexpressed as 
polynomials in R^ m . For even N, the irreducible zero modes with the lowest scaling 
dimension have the form 

f„ a (R) = RlRl .R 2 + [...] (81) 

J N,0\ / 12 34 (JV-l)JV 1 L J V J 

where [ . . . ] denotes a combination of terms that depend on the positions of (N — 1) or less 
particles. For four particles, the zero mode is 2(d + 2)i? 1 2 2 i? 3 2 4 — d(R* 2 + R34), our recurrent 
example. The terms [ . . . ] are not uniquely determined since any degree N zero mode for a 
smaller number of points might be added. Furthermore, permutations of the points in f N0 
give other zero modes so that we may symmetrize the above expressions and look only at 
the permutation-invariant modes. The scaling dimension ( N0 of f N0 is clearly equal to N. 
This linear growth signals the absence of attractive effects between the particles diffusing 
with a constant diffusivity (no particle binding in the shape evolution). As we shall see 
in Sect. III. C.l, this leads to the disappearance of the intermittency in the advected scalar 
field, that becomes Gaussian in the limit £ — > 0. 

To the linear order in £, the operator M. N will differ from the Laplacian by a second 
order differential operator — £ V, involving logarithmic terms oc ln(r nm ). The zero mode and 
its scaling dimension are expanded as fo + £ f\ and N + ££1, respectively. The lowest order 
term / is given by the symmetrization of (81). As usual in such problems, the degeneracy 
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hidden in [ . . . ] may be lifted by the perturbation that fixes fo for each zero mode, see 
below. At the first order in £, the equations that define the zero modes and their scaling 
dimension reduce to the relations 



V 2 A = Vf , (Rd R ~ N) fi = Clio- (82) 

Given an arbitrary zero mode / , one shows that the first equation admits a solution of the 
form fi = h+ h nm ln(r nm ) with 0(<i)-invariant, degree N polynomials h nm and h, the 
latter being determined up to zero modes of V 2 . Note that the function (Rd R — N) fi — 
h nm is also annihilated by the Laplacian: 



V 2 (it> ,9, - AO A = (Rd R - N + 2)V 2 h = (Rd R - N + 2)Vf 

= ([Rd R ,V]+2V)f + V(Rd R -N)f = 0. (83) 

The last equality follows from the scaling of / and the fact that the commutators of R d R 
with M N and V are (£ — 2) M N and — V 2 — 2 V, respectively. One obtains this way a linear 
map T on the space of the degree N zero modes of the Laplacian: T f = (Rd R — N) f\. 
The second equation in (82) states that f must be chosen as an eigenstate of the map T. 
Furthermore, the function should not belong to the subspace of unit codimension of the zero 
modes that do not depend on all the points. It is easy to see that such subspace is preserved 
by the map T. As the result, the eigenvalue Ci is equal to the ratio between the coefficients 
of _R 2 2 _R 2 4 . . . R^ N _ 1)N in r/ and in f . The latter is easy to extract, see (Bernard et al, 
1996) for the details, and yields the result £i = — N 2 ^+i) or > equivalently, 

C, - f(2-« - yg^t + 0<a (84, 

giving the leading correction to the scaling dimension of the lowest irreducible zero mode. 
Note that to that order ( N0 is a concave function of N. Higher-order terms in £ have been 
analyzed in (Adzhemyan et al, 1998) (the second order) and in (Adzhemyan et al, 2001) 
(the third order). The latter papers used a renormalization group resummation of the small £ 
perturbative series for the correlation functions of the scalar gradients in conjunction with an 
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operator product expansion, see Sec. III.C.l. The expression (84) may be easily generalized 
to the compressible Kraichnan ensemble of compressibility degree p. The correction d for 
the tracer exponent picks up an additional factor (l + 2p) (Gawedzki and Vergassola, 2000). 
Higher order corrections may be found in (Antonov and Honkonen, 2000). The behavior of 
the density correlation functions was analyzed in (Adzhemyan and Antonov, 1998; Gaw§dzki 
and Vergassola, 2000; Antonov and Honkonen, 2000). 

For large dimensionality d, it is convenient to use the variables x nm = Rnm as the 
independent coordinates 6 to make the (i-dependence in A4 N explicit. Up to higher orders 
in 1/d, the operator M N oc C— \ U, where £ = d~ l [(d— l)d x + (2 — £)x nm d 2 } and 

n<m nm Xnm 

U is a second order (/-independent differential operator mixing derivatives over different 
x nm . We shall treat C as the unperturbed operator and —\U as a perturbation. The 
inclusion into C of the diagonal terms oc | makes the unperturbed operator of the same 
(second) order in derivatives as the perturbation and renders the perturbative expansion 
less singular. The irreducible zero modes of C with the lowest dimension are given by an 
expression similar to (81): 

f Nfi {R) = x 12 x. i4 ... x (N ^ 1)N + [...], (85) 

and the permutations thereof. Their scaling dimension is y(2 — £). For N = 4, one may 
for example take f^o = x 12 x 34 — 2 (2-t) ( x i2 + X U' ^ s m ^ ne £- ex P ans i° n ; i n order to take 
into account the perturbation U, one has to solve the equations 

Ch = Uf , (E^-t)/! = ^/o- (86) 

n<m 

One checks again that r/ = ( x nm d x — y) fi is annihilated by C. In order to calculate 
(i, it remains to find the coefficient of x 12 . . -x, N _ 1)N in T/o- In its dependence on the x's, 
the function Ufo scales with power (y — 1). One finds fi by applying the inverse of the 
operator C to it. r/ is then obtained by gathering the coefficients of the logarithmic terms 

6 Their values are restricted only by the triangle inequalities between the interparticle distances. 
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in fi, see (Chertkov et al, 1995b) for the details. When d — > oo, the operator C reduces 
to the first order one CJ = Y) d . This signals that the particle evolution becomes 

n<m nm 

deterministic at d = oo, with all x nm growing linearly in time. If one is interested only in 
the ^ correction to the scaling exponent and not in the zero mode, then it is possible to 
use directly the more natural (but more singular) decomposition Ai N oc £ — \ U' . The 
leading zero modes of £ also have the form (85). Noting that £ is a translation operator, 
the zero mode r'/ may be obtained as the coefficient of the logarithmically divergent term 

(•OO 

in / U'fo(x nm — t) dt , see (Chertkov and Falkovich, 1996). In both approaches, the final 
Jo 

result is 

C^fp-O-^^f + O^), (87) 

which is consistent with the small £ expression (84). 

The non-isotropic zero modes, as well as those for odd N, may be studied similarly. The 
zero modes of fixed scaling dimension form a representation of the rotation group SO(d) 
which may be decomposed into irreducible components. In particular, one may consider the 
components corresponding to the symmetric tensor products of the defining representation 
of SO(d) , labeled by the angular momentum j (the multiplicity of the tensor product). 
For two particles, no other representations of SO(d) appear. The 2-point operator J\4 2 
becomes in each angular momentum sector an explicit second-order differential operator in 
the radial variable. It is then straightforward to extract the scaling dimensions of its zero 
modes: 

Note that ^ = 1 in any d corresponding to linear zero modes. For the 3-point operator, 
the lowest scaling dimensions are (® Q = 4 — ^5p£ + 0(£ 2 ) (Gat et al, 1997a,b) and 
Ca.o = 3 - f + 0(e) (Pumir, 1996 and 1997) or = 3 - £ - f + 0{±) (Gutman 
and Balkovsky, 1996). For even iV and j, Arad, L'vov et al. (2000) have obtained the 
generalization of (84) in the form 
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a N (9 n f N ( N ~ 2 ) j(j + d-2)(d + i) \ 2 

C - = T (2 "° " {-2(dTW " 2(d + 2)(d-l) + ^ } - (89) 

The effective expansion parameter in the small £ or large d approach turns out to be 
^ ^ d so that neither of them is applicable to the region of the almost smooth velocity fields. 
This region requires a different perturbative technique exploiting the numerous symmetries 
exhibited by the multiparticle evolution in the limiting case £ = 2. Those symmetries were 
first noticed and employed to derive an exact solution for the zero modes by Shraiman and 
Siggia (1996). The expression of the multiparticle operators at £ = 2 reads 

M N = D 1 (dH 2 -(d+l)J 2 ) , (90) 

with H 2 = J2 HijHji and J 2 = — J2 Jjj denoting the Casimir operators of the group 

ij i<j 

SL(d) and of its SO(d) subgroup acting on the index i — 1, . . . , d of the particle positions 
r l n . The corresponding generators are given by — ( — r% nS r i + ^«i( r n^ an d 
Jij = Hij — Hji. The relation (90), that may be easily checked directly, is consistent with 
the expression (68) for the heat kernel of M N . As mentioned in Sect. II.B(i), the right 
hand side of (90) is indeed the generator of the diffusion process W(t) on the group SL(d). 
In their analysis, Shraiman and Siggia (1995) employed an alternative expression for the 
multiparticle operators, exhibiting yet another symmetry of the smooth case: 

M N = D 1 (dG 2 -(d+l)J 2 + d -^A(A + d N )^ , (91) 

N-l 

where G 2 = G nm G mn is the quadratic Casimir of SL(N — 1) acting on the index 

n,rn=l 

n = 1, . . . , N — 1 of the difference variables r % and A = Y\ r % n V r i is the generator of 

i, n 

the overall dilations. For three points, one may then decompose the scaling translationally 
invariant functions into the eigenfunctions of G 2 , L 2 and other generators commuting with 
the latter and with A. The zero modes of A4 3 at £ = 2 have the lowest scaling dimension 
equal to unity and vanishing in the angular momentum sectors j — 1 and j = , see (Pumir 
et al, 1997) and (Balkovsky et al, 1997a). They have infinite multiplicity since their space 
carries an infinite-dimensional representation of SL{2). Similar zero modes exist for any 
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higher scaling dimension. In 2d, for example, the scaling dimension of a three-particle zero 
mode may be raised multiplying it by a power of det(r* 3 ). The continuous spectrum of the 
dimensions and the infinite degeneracy of the zero modes in the smooth case is one source 
of the difficulties. Another related difficulty is that for £ = 2 the principal symbol of M3 
looses positive-definiteness not only when two of the points coincide, but also when all the 
three points become collinear. That leads to the domination for quasi-collinear geometries of 
the perturbative terms in M. z over the unperturbed ones. The problem requires a boundary 
layer approach developed first by Pumir et al. (1997) for the j — 1 sector with the conclusion 
that the minimal scaling dimension of the zero mode behaves as 1 + o(2 — £). Similar 
techniques led Balkovsky et al. (1997a) to argue for a 0{\J2 — £) behavior of the minimal 
scaling dimension of the isotropic zero modes. The three-particle zero mode equation was 
solved numerically for the whole range of values of £ by Pumir (1997) for j = 1 and d = 2, 3 
and by Gat et al. (1997a,b) in the isotropic case for d — 2, 3, 4 . Their results are compatible 
with the perturbative analysis around £ = and £ = 2 , with a smooth interpolation for the 
intermediate values (no crossing between different branches of the zero modes). Analytical 
non-perturbative calculations of the zero modes were performed for the passive scalar shell 
models, where the degrees of freedom are discrete. We refer the interested reader to the 
original works (Benzi et al, 1997; Andersen and Muratore-Ginanneschi, 1999) and to (Bohr 
et al, 1998) for an introduction to shell models. 



III. PASSIVE FIELDS 

The results on the statistics of Lagrangian trajectories presented in Chapter II will be 
used here to analyze the properties of passively advected scalar and vector fields. The 
qualification "passive" means that we disregard the back-reaction of the advected fields on 
the advecting velocity. We shall treat both a scalar per unit mass (a tracer field), satisfying 
the equation 

d t 9 + v ■ V9 = kV 2 6 , (92) 
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and the density per unit volume, whose evolution is governed by 

d t n + V • (n v) = nV 2 n . (93) 

For incompressible flows, the two equations are obviously coinciding. Examples of passive 
vector fields are provided by the gradient of a tracer u> = V0, obeying 

d t u + V(u • w) = kV 2 w, (94) 

and the divergenceless magnetic field evolving in incompressible flow according to 

d t B + v VB — B V« = kV 2 B . (95) 

Two broad situations will be distinguished : forced and unforced. The evolution equations 
for the latter are (92)-(95). They will be analyzed in Sect. III. A. For the former, a pumping 
mechanism such as a forcing term is present and steady states might be established. The 
rest of the Chapter treats steady cascades of passive fields under the action of a permanent 
pumping. As we shall see below, the advection equations may be easily solved in terms of 
the Lagrangian flow, hence the relation between the behavior of the advected fields and of 
the fluid particles. In particular, the multipoint statistics of the advected fields will appear 
to be closely linked to the collective behavior of the separating Lagrangian particles. An 
introduction to the passive advection problem may be found in (Shraiman and Siggia, 2000). 

A. Unforced evolution of passive scalar and vector fields 

The physical situation of interest is that the initial passive field or its distribution is 
prescribed and the problem is to determine the field distribution at a later time t. The 
simplest question to address is which fields have their amplitudes decaying in time and which 
growing, assuming the velocity field to be statistically steady. A tracer field always decays 
because of dissipative effects, with the law of decay depending on the velocity properties. The 
fluctuations of a passive density may grow in a compressible flow, with this growth saturated 
by diffusion after some time. The fluctuations of both u> and B may grow exponentially 
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as long as diffusion is unimportant. After diffusion comes into play, their destinies are 
different : u> decays, while the magnetic field continues to grow. This growth is known 
as dynamo process and it continues until saturated by the back-reaction of the magnetic 
field on the velocity. Another important issue here is the presence or absence of a dynamic 
self-similarity : for example, is it possible to present the time-dependent PDF V(9; t) as a 
function of a single argument? In other words, does the form of the PDF remain invariant in 
time apart from a rescaling of the field? We shall show that for large times the scalar PDF 
tends to a self-similar limit when the advecting velocity is nonsmooth, while self-similarity 
is broken in smooth velocities. 

1. Backward and forward in time Lagrangian description 

If the advecting velocities are smooth and if the diffusive terms are negligible, 7 the 
advection equations may be easily solved in terms of the Lagrangian flow. To calculate 
the value of a passively advected field at a given time one has to trace the field evolution 
backwards along the Lagrangian trajectories. This is to be contrasted with the description of 
the particles in the previous Chapter, which was developed in terms of the forward evolution. 
The tracer 9 stays constant along the Lagrangian trajectories: 

9(r,t) = 9(R(0;r,t),0), (96) 

where R( ■ ;r,t) denotes the Lagrangian trajectory passing at time t through the point r. 
The density n changes along the trajectory as the inverse of the volume contraction factor. 
Let us consider the matrix W(t; r) = W(t; R(0; r, t)) , where W(t; r) is given by (13) and 
describes the forward evolution of small separations of the Lagrangian trajectories starting 
at time zero near r. The volume contraction factor is det(W(t; r)) and 

7 Recall that the Schmidt number u/k, also called Prandtl number when considering temperature 
or magnetic fields, is assumed to be large. 
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n(r,t) = [detiWfcr))]- 1 n(R(0;r,t),0) . (97) 

Note that the matrix W(t; r) is the inverse of the backward-in-time evolution matrix 
W'(t; r) with the matrix elements dK l (0;r,t)/dr j . This is indeed implied by the identity 
R(t; R(0; r, t), 0) = r and the chain rule for differentiation. The solution of the evolution 
equation for the gradient of the tracer is obtained by differentiating (96): 

w(r,t) = (W(*;r)-T"WW),0). (98) 

Finally, the magnetic field satisfies 

B(r,t) = W(t;r) B(R(0;r,t),0) . (99) 

The relations (96) to (99) give, in the absence of forcing and diffusion, the solutions of the 
initial value problem for the advection equations in a given realization of a smooth velocity. 

For non-zero k, the solutions of the scalar equations are given essentially by the same 
expressions. However, R( • ;r,t) denotes now a noisy Lagrangian trajectory satisfying the 
stochastic equation (5) and passing through r at time t and the right hand sides of the 
equations (96) to (99) should be averaged over the noise using the Ito formula of stochastic 
calculus discussed in Appendix. These solutions may be rewritten using the transition PDF's 
p(r, s; R, 1 1 v) introduced in Sect. II. C, see (43) and describing the probability density to 
find the noisy particle at time t at position R, given its time s position r. One has 

0(r, t)= J p(r, t;R,0\ v) 0(R, 0) dR , n(r, t) = J p{R, 0;r,t\ v) n(R, 0) dR . (100) 

The two PDF's appearing in these formulae, one backward and the other forward in time, 
coincide for incompressible velocities but they are generally unequal for the compressible 
cases. For nonsmooth velocities, those PDF's continue to make sense and we shall use (100) 
to define the solutions of the scalar advection equations in that case. As for the vector 
fields, their properties depend both on the noisy Lagrangian trajectory endpoints and on 
the matrices W(t; r), that are well defined only in smooth velocities. The formal procedure 
for nonsmooth velocities is to first impose a viscous cutoff, smoothing the velocity behavior 
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at small scales, and then removing it. When this is done, some properties of the field remain 
well defined and may be analyzed (see, for instance, Section III.C.3). 

In random velocity fields, the advected quantities become random fields whose statistics 
may be probed by considering the equal-time correlation functions. In particular, those of 
a tracer evolve according to 

C N (r,t) = (9( ri ,t) ...9(r N ,t)) = J V N (r ; R\-t) 0(^,0) ... 9(R N ,0) dR. (101) 

Here, as in Sect. II. C, the Green functions V N are the joint PDF's of the equal-time po- 
sitions of iV fluid particles, see (65). For the correlators of the density n, the backward 
propagator in (101) should be replaced by its forward version, in agreement with (100). If 
the initial data are random and independent of the velocities, they may be easily averaged 
over. For a tracer with a Gaussian, mean zero initial distribution: 

C 2n (r, t) = J V 2n (r; R, -t) [ C 2 (R 12 , 0) . . . C 2 (i* (2n _i) 2n , o) + . . . ] dR , (102) 

where, according to the Wick rule, the dots stand for the other pairings of the 2n points. 

Let us now briefly discuss the compressible case, where the statistics of the matrices 
W and W generally do not coincide. As we have already discussed in Section II. D, every 
trajectory then comes with its own weight determined by the local rate of volume change 
and exhibited by the Lagrangian average of a function f(W): 

J f(W(t;r)) % = J f(W(t;R)) det(W(t;R)) f . (103) 

The relation in (103) simply follows from the definition of W. The volume change factor 
det(VF) = expX^Pi, with the same notation as in Sect. II. B. Recall that only the average 
of the determinant is generally equal to unity for compressible flow. The averages of the 
S'0((i)-invariant functions of W are described for large times by the large-deviation function 
H = H — Pi/t, with the last term coming from the volume factor. The corresponding 
Lyapunov exponents \ are determined by the extremum of H (Balkovsky et al, 1999a). 
The exponents generally depend on the form of the entropy function H and cannot be 
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expressed via the Lyapunov exponents Aj only. Since the matrix W of the backward 
evolution is the inverse of W, the backward Lyapunov exponents are given by — \d-i+i and 
not by the naive guess —\d-i+i- In particular, the interparticle distance diverges backward in 
time with the exponent — \d- The same way as we have shown in Sect. II. D that J2 A« < in 
a compressible flow, one shows that X) Aj > (implying Ai > 0). For the forward Lagrangian 
evolution we thus have an average compression of volumes, whereas passive fields rather feel 
an average expansion. Indeed, as we go away from the moment where we imposed a uniform 
Lagrangian measure, the rate of change of the volume is becoming negative in a fluctuating 
compressible flow. 

The forward and backward in time Lyapunov exponents coincide if the statistics of the 
velocities is time-reversible, i.e. if v{r,t) and —v{r,—t) are identically distributed. More 
generally, the entire distributions of the forward and the backward in time stretching rates 
coincide in that case: 

H( Pl /t - Ai, . . . , Pd /t - X d ) = H(- Pd /t - Ai, . . . , -pi/f - A d ) + 2 Pi/t • (104) 

i 

This is an example of the time-reversibility symmetry of the large deviation entropy function 
that was described by Gallavotti and Cohen (1995). The symmetry holds also for a 5- 
correlated strain, although the above finite-volume argument (103) does not apply directly 
to this case. Recall that in that instance the entropy function (59) describes the large 
deviations of the stretching rates of the matrix W(t) given by the Ito version of (16). For 
the inverse evolution, the strain a(s) should be replaced by cr'(s) = —a(t — s) and the 
matrix W'(t) has the same distribution as W(t). The matrix W(t) = W'^)^ 1 is then 
given by (16) with the anti-Ito regularization and the relation between the conventions (see 
A4) implies W(t) = W(t) e~ 2 S Xit / d . Realistic turbulent flows are irreversible because of 
the dissipation so that the symmetry (104), that was confirmed in an experimental situation 
(Ciliberto and Laroche, 1998), may be at most approximate. 
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2. Quasi-Lagrangian description of the advection 

Important insights into the advection mechanisms are obtained by eliminating global 
sweeping effects and describing the advected fields in a frame whose origin moves with 
the fluid. This picture of the hydrodynamic evolution, known under the name of quasi- 
Lagrangian description, is intermediate between the static Eulerian and the dynamic La- 
grangian ones (Monin, 1959; Belinicher and L'vov, 1987). Specifically, quasi-Lagrangian 
fields are defined as 



where ip stands for any Eulerian field, scalar or vector, introduced previously and R(t; r , 0) 
is the Lagrangian trajectory passing through r at time zero. The quasi-Lagrangian fields 
satisfy the same evolution equations as the Eulerian ones except for the replacement of the 
advective term by [v(r, t) — v(0, t)] • V. If the incompressible velocity and the initial values 
of the advected field are statistically homogeneous, the equal-time statistics of the quasi- 
Lagrangian and the Eulerian fields coincide. The equal-time statistics is indeed independent 
of the initial point r . The equality of the equal-time Eulerian and quasi-Lagrangian distri- 
butions follows then by averaging first over r and then over the velocity and by changing the 
variables r i— > R(t; r , 0) . The equality does not hold if the initial values of the advected 
fields are non- homogeneous. 

It will be specially convenient to use the quasi-Lagrangian picture for distances r much 
smaller than the viscous scale, i.e. in the Batchelor regime (Falkovich and Lebedev, 1994). 
The variations in the velocity gradients may then be ignored so that v(r, 0) — v(0, t) ps a(t)r. 
In this case, the velocity field enters into the advection equations only through the time 
dependent strain matrix. For the tracer, one obtains then the evolution equation 



This may be solved as before using now the noisy Lagrangian trajectories for a velocity 
linear in the spatial variables: 



4>(r,t) 



i>(r + R(t; r o ,0),f), 



(105) 




(106) 
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§(r,t) = 0(w(t)- l r -V2^ ^V(s)^ 1 rf/3(s), 0), 



(107) 



where W(t) is the evolution matrix of (16). The overline denotes the average over the noise 
which is easily performed for incompressible velocity fields using the Fourier representation: 



0(r, t) = J 9{W{t) T k, 0) exp [i k ■ r - k ■ Q(t) k 

with Q(t) = k t Wit^MisfWis^Witf . 
Jo 



dk 

(2w) d 



(108) 
(109) 



3. Decay of tracer fluctuations 

For practical applications, e.g. in the diffusion of pollution, the most relevant quantity 
is the average {9(r, t)}. It follows from (101) that the average concentration is related to the 
single particle propagation discussed in Sect. II. A. For times longer than the Lagrangian 
correlation time, the particle diffuses and (9) obeys the effective heat equation 

d t (6(r,t))=D^V i V j (e(r,t)), (110) 

with the effective diffusivity D' l J given by (9). The decay of higher-order moments and 
multipoint correlation functions involves multiparticle propagation and it is sensitive to the 
degree of smoothness of the velocity field. 

The simplest decay problem is that of a uniform scalar spot of size £ released in the fluid. 
Another relevant situation is that where a homogeneous statistics with correlations decay- 
ing on the scale £ is initially prescribed. The corresponding decay problems are discussed 
hereafter for the two cases of smooth and nonsmooth incompressible flow, 
i) Smooth velocity. Let us consider an initial scalar configuration given in the form 
of a single spot of size £. Its average spatial distribution at later times is given by the 
solution of (110) with the appropriate initial condition. On the other hand, the decay of the 
scalar in the spot as it is carried with the flow corresponds to the evolution of 9 , defined 
by (105). We assume the Schmidt/Prandtl number v/k to be large so that the viscous 
scale i] is much larger than the diffusive scale r^. We shall be considering the 3d situation, 
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where the diffusive scale = ^ and A 3 is the most negative Lyapunov exponent 
defined in Section II.B.f . As shown there, for times t <C £<f — | A3 1 1 ln(£/r^), the diffusion is 
unimportant and the values of the scalar field inside the spot do not change. At later times, 
the width in the direction of the negative Lyapunov exponent A3 is frozen at r^, while the 
spot keeps growing exponentially in the other two directions 8 . The freezing of the contracting 
direction at r d thus results in an exponential growth of the volume oc exp(p x + p 2 ). Hence, 
the scalar moment of order a measured at locations inside the spot will decrease as the 
average of exp[— a(p 1 + p 2 )]. The resulting decay laws exp(— j a t) may be calculated using 
the PDF (27) of the stretching variables pj. More formally, the scalar moments inside the 
spot are captured by the quasi-Lagrangian single-point statistics. Following Balkovsky and 
Fouxon (1999), let us take in (108) a Gaussian initial configuration 9(k,0) = exp[— \{ik) 2 }. 
As a result, 

0(0, f) = y exp [- \l 2 k-I(t)k] ^3 oc det J(t)~ 1/2 = e~^ p \ (111) 

where I(t) is the mean tensor of inertia introduced in Sect. II.B.l, see (25). Using the PDF 
(27), one obtains then 

(9 a (t)) oc J exp[-a( Pl + p 2 )-tH( Pl /t- \ 1 ,p 2 /t- \ 2 )] d Pl dp 2 . (112) 

At large times, the integral is determined by the saddle point. At small a, it lies within 
the parabolic domain of H and the decay rate 7 a increases quadratically with the order a. 
At large enough orders, the integral is dominated by the rare realizations where the volume 
of the spot does not grow in time and the growth rates become independent of the order 
(Shraiman and Siggia, 1994; Son, 1999; Balkovsky and Fouxon, 1999). That conclusion is 
confirmed experimentally (Groisman and Steinberg, 2000). 

An alternative way to describe the decay of (9 N (t)) is to take iV fluid particles that come 
to the given point r at time t and to track them back to the initial time. The realizations 



8 As in Section II. B, we consider the case of two non-negative Lyapunov exponents; the arguments 
are easily modified for two non-positive exponents. 
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contributing to the moments are those for which all the particles were initially inside the 
original spot of size £, see (107). Looking backward in time, we see that the molecular 
noise splits the particles by small separations of order 0(rd) during a time interval of the 
order td = r\j k (A^ 1 ) near t. After that, the advection takes over. The realizations 
contributing to the moments are those with the interparticle separations almost orthogonal 
to the (backward) expanding direction p 3 of W~ x . More exactly, they should form an angle 
~ (^/Vd)e P3 with the plane orthogonal to the expanding direction. Such separations occupy 
a solid angle fraction of the same order. Since we now track particles moving due to the 
advection (the molecular noise is accounted for by the finite splitting) then J2 Pi = and 
(112) follows. 

The same simple arguments lead to the result for the case of a random initial condition 
with zero mean. Let us first consider the case when it is Gaussian with correlation length £. 
It follows from (102) that the realizations contributing to (9 2n ) = (0 2 ") are those where n 
independent pairs of particles are separated by distances smaller than £ at time t — 0. The 
moments are therefore given by 

(6 2n (t)) oc jexp[-n(p 1 +p 2 )-tH(p 1 /t-X 1 ,p 2 /t-X 2 )] d Pl dp 2 . (113) 

Note that the result is in fact independent of the scalar field initial statistics. Indeed, for a 
non-Gaussian field we should average over the Lagrangian trajectories the initial correlation 
function C 2n (R(0),0) that involves a non-connected and a connected part. The latter is 
assumed to be integrable with respect to the 2n — 1 separation vectors (and thus to depend 
on them). Each dependence brings an exp [—pi — p 2 ] factor and the connected part will thus 
give a subleading contribution with respect to the non-connected one. The above results 
were first obtained by Balkovsky and Fouxon (1999) using different arguments (see Sections 
III. 4 and III. 5). Remark the square root of the volume factor appearing in (113) as distinct 
from (112). In the language of spots, this is explained by the mutual cancellations of the 
tracer from different spots and the ensuing law of large numbers. Indeed, different blobs of 
size £ with initially uncorrelated values of the scalar will overlap at time t and the r.m.s. 
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value of 9 will be proportional to the square root of the number of spots oc exp(pi + p 2 ). The 
same qualitative conclusions drawn previously about the decay rates 7 Q may be obtained 
from (113). In particular, Balkovsky and Fouxon (1999) performed the explicit calculation 
for the short-correlated case (28). The result is 

7o = - Urn il„(|r) cx o(l- |), (IM) 

for a < 4 and 7 a = const, for a > 4. An important remark is that the PDF of the decaying 
scalar is not self-similar in a smooth velocity field. The PDF is indeed becoming more and 
more intermittent with time, as signaled by the growth of the kurtosis (\6\ a )/(6 2 ) ^ for 



a > 2. The previous arguments may be easily generalized to the case of compressible flow, 
ii) Nonsmooth velocity. For the decay in incompressible nonsmooth flow, we shall 
specifically consider the case of a time reversible Kraichnan velocity field. The comments on 
the general case are reserved to the end of the section. The simplest objects to investigate 
are the single-point moments (6 2n {t)^ and we are interested in their long-time behavior 
t 3> l 2 ~^/Di. Here, £ is the correlation length of the random initial field and D 1 enters the 
velocity 2-point function as in (48). Using (101) and the scaling property (71) of the Green 
function we obtain 

(e 2n (t)) = J V 2n (0; R, -1) C 2n (t^R, 0) dR. (115) 

There are two universality classes for this problem, corresponding to either non-zero or 
vanishing value of the so-called Corrsin integral J = / C 2 {r,t) dr. Note that the integral is 
generally preserved in time by the passive scalar dynamics. 

We concentrate here on the case J 7^ and refer the interested reader to the original 

d 1 

paper by Chaves et al. (2001) for more details. For J 7^ 0, the function t 2 -? C*2(t 2 -«r,0) 
tends to Jo5(r) in the long-time limit and (115) is reduced to 

(9 2n (t)) » (2n - 1)!! f V 2n (0; R u R u R n , R n , -1) dR, (116) 

for a Gaussian initial condition. A few remarks are in order. First, the previous formula 
shows that the behavior in time is self-similar. In other words, the single point PDF V(t, 6) 
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takes the form t^^Q^t^^O). That means that the PDF of Qj\fl is asymptotically time- 
independent as was hypothetized by Sinai and Yakhot (1989), with e(t) = k,((V6) 2 ) being 
time-dependent (decreasing) dissipation rate. This should be contrasted with the lack of 
self-similarity found previously for the smooth case. Second, the result is asymptotically 
independent of the initial statistics (of course, within the universality class Jo 7^ 0). As in 
the previous subsection, this follows from the fact that the connected non-Gaussian part 
of C 2n depends on more than n separation vectors. Its contribution is therefore decaying 

nd 

faster than t 2-5 . Third, it follows from (116) that the long-time PDF, although universal, 
is generally non-Gaussian. Its Gaussianity would indeed imply the factorization of the 
probability for the 2n particles to collapse in pairs at unit time. Due to the correlations 
existing among the particle trajectories, this is generally not the case, except for £ = where 
the particles are independent. The degree of non- Gaussianity is thus expected to increase 
with £, as confirmed by the numerical simulations presented in (Chaves et al, 2001). 

Other statistical quantities of interest are the structure functions S 2n (r,t) = ([6(r,t) — 
9(0,t)] 2n ) related to the correlation functions by 

S 2n (r,t) = j\..j o 1 d fli ...d, 2n C 2n (fi 1 r,..., hn r,t) JJd^ = A(r) C 2n {-) . (117) 

To analyze their long-time behavior, we proceed similarly as in (115) and use the asymptotic 
expansion (75) to obtain 



S 2n (r,t) = j A(r 3=?r) V 2n ( • ; R; -1) C 2n {t~ R,0) dR 

- A{r)f2n ' oi ' ) jg 2 n,o(R,-^C 2n (t^R,o)dR oc f^V'Vm (118) 



C2n,0 



Here, f 2n ,o is the irreducible zero mode in (75) with the lowest dimension and the scalar 
integral scale £(t) oc t~ . As we shall see in Sect. III.C.l, the quantities ( 2n = ( 2ri:0 give 
also the scaling exponents of the structure function in the stationary state established in the 
forced case. 

Let us conclude this subsection by briefly discussing the scalar decay for velocity fields 
having finite correlation times. The key ingredient for the self-similarity of the scalar PDF 

72 



is the rescaling (71) of the propagator. Such property is generally expected to hold (at least 
for large enough times) for self-similar velocity fields regardless of their correlation times. 
This has been confirmed by the numerical simulations in (Chaves et al, 2001). For an 
intermittent velocity field the presence of various scaling exponents makes it unlikely that 
the propagator possesses a rescaling property like (71). The self-similarity in time of the 
scalar distribution might then be broken. 

4. Growth of density fluctuations in compressible flow 

The evolution of a passive density field n(r, t) is governed by the equation (93). In smooth 
velocities and in the absence of diffusion, its solution is read from (97), where we shall take 
the initial field on the right hand side to be uniform. This gives n(r,t) = [det(W(t; r))]^ 1 . 
Performing the velocity average and recalling the long-time asymptotics of the W statistics, 
we obtain 



The moments at long times may be calculated by the saddle-point method and they are 
generally behaving as oc exp(7 a t). The growth rate function 7 a is convex, due to Holder 
inequality, and vanishes both at the origin and for a — 1 (by the total mass conservation). 
This leads to the conclusion that 7 a is negative for < a < 1 and is otherwise positive: 
low-order moments decay, whereas high-order and negative moments grow. For a Kraichnan 
velocity field, the large deviations function H is given by (59) and the density field becomes 
lognormal with j a oc a(a — 1) (Klyatskin and Gurarie, 1999). Note that the asymptotic 
rate (\nn)/t is given by the derivative at the origin of j a and it is equal to — J2 Aj < 0. 
The density is thus decaying in almost any realization if the sum of the Lyapunov exponents 
is nonzero. Since the mean density is conserved, it has to grow in some (smaller and 
smaller) regions, which implies the growth of high moments. The amplification of negative 
moments is due to the growth of low density regions. The positive quantity — J2 \ has the 
interpretation of the mean (Gibbs) entropy production rate per unit volume. Indeed, if we 




(1 - a)Y,Pi ~tH(p 1 /t - Ai, . . .,p d /t - X d ) Y[dpi . 



(119) 
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define the Gibbs entropy S(n) as — f(kin)ndr then, by (103), S(n) = / In det(W(t; r)) dr. 
Since In det(W) = Pi, the entropy transfered to the environment per unit time and unit 
volume is — J2Pi/t an d it is asymptotically equal to — J2^i > 0, see (Ruelle, 1997). 

The behavior of the density moments discussed above is the effect of a linear but random 
hyperbolic stretching and contracting evolution (13) of the trajectory separations. In a finite 
volume, the linear evolution is eventually superposed with non-linear bending and folding 
effects. In order to capture the combined impact of the linear and the non-linear dynamics 
at long times, one may observe at fixed time t the density produced from an initially uniform 
distribution imposed at much earlier times to- When to — > — oo and if Ai > 0, the density 
approaches weakly, i.e. in integrals against test functions, a realization-dependent fractal 
density n*(r, t) in almost all the realizations of the velocity. The resulting density field is the 
so-called SRB (Sinai- Ruelle-Bowen) measure, see e.g. (Kiffer, 1988). The fractal dimension 
of the SRB measures may be read from the values of the Lyapunov exponents (Frederikson 
et ai, 1983). For the Kraichnan ensemble of smooth velocities, the SRB measures were 
first discussed by (Le Jan, 1985). In 2d, they have a fractal dimension equal to 1 + if 
< p < \ . In 3d, the dimension is 2 + ±=§| if < p < § and 1 + ^2 if § < p < §, 
where p is the compressibility degree. 

The above considerations show that, as long as one can neglect diffusion, the passive 
density fluctuations grow in a random compressible flow. One particular case of the above 
phenomena is the clustering of inertial particles in an incompressible turbulent flow, see 
(Balkovsky et ai, 2001) where the theory for a general flow and the account of the diffusion 
effects that eventually stops the density growth were presented. 

5. Gradients of the passive scalar in a smooth velocity 

For the passive scalar gradients uj = V# in an unforced incompressible situation, the 
equation to be solved is (94). The initial distribution is assumed statistically homogeneous 
with a finite correlation length. As discussed previously, one may treat diffusion either by 
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adding a Brownian motion to the backward Lagrangian trajectories or by using the Fourier 
transform method (108). For pedagogical reasons, we choose here the latter and solve 
(94) by simply taking the gradient of the scalar expression (108). The long-time limit is 
independent of the initial scalar statistics (Balkovsky and Fouxon, 1999) and it is convenient 
to take it Gaussian with the 2-point function oc exp[— ^(r/f) 2 ]. The averaging over the 
initial statistics for the generating function Z{y) = (exp [iy ■ u>]) reduces then to Gaussian 
integrals involving the matrix I(t) determined by (25). The inverse Fourier transform is 
given by another Gaussian integral over y and one finally obtains for the PDF of u> : 



As may be seen from (25), during the initial period t < t d = |A^ X | \n(£/r d ), the diffusion is 
unimportant, the contribution of the matrix Q to / is negligible, the determinant of the 
latter is unity and uo 2 grows as the trace of J -1 . In other words, the statistics of lnu; and 
of — p d coincide in the absence of diffusion. The statistics of the gradients can therefore 
be immediately taken over from Section II. B. The growth rate (2t) _1 (lnu; 2 ) approaches 
| Ad | while the gradient PDF depends on the entropy function. For the Kraichnan model 
(28), the PDF is lognormal with the average Did(d — l)t and the variance 2Di(d — l)t 
read directly from (32). This result was obtained by Kraichnan (1974) using the fact that, 
without diffusion, u satisfies the same equation as the distance between two particles, whose 
PDF is given by (23). 

As time increases, the wavenumbers (evolving as k = a T k ) reach the diffusive scale r^ 1 
and the diffusive effects start to modify the PDF, propagating to lower and lower values of 
uo. High moments first and then lower ones will start to decrease. The law of decay at t 3> t d 
can be deduced from (120). Considering this expression in the eigenbasis of the matrix /, 
we observe that the dominant component of uj coincides with the largest eigendirection of 
the J -1 matrix, i.e. the one along the p d axis. Recalling from the Section II. B that the 
distribution of p d is stationary, we infer that ^|cj| a (t)^ oc ^(det/)~ a / 4 ^} . The comparison 
with (113) shows that the decay laws for the scalar and its gradients coincide (Son, 1999; 




(120) 
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Balkovsky and Fouxon, 1999). This is qualitatively understood by estimating cu ~ 9/£ min , 
where £ min is the smallest size of the spot. Noting that 9 and £ min are independent and that 
£ m i n ~ e pd £ at large times has a stationary statistics concentrated around r<j, it is quite clear 
that the decrease of u is due to the decrease of 9. 

6. Magnetic dynamo 

In this Section, we consider the generation of inhomogeneous magnetic fluctuations below 
the viscous scale of incompressible turbulence. The question is relevant for astrophysical 
applications as the magnetic fields of stars and galaxies are thought to have their origin in the 
turbulent dynamo action (Moffatt, 1978; Parker, 1979; Zeldovich et al, 1983; Childress and 
Gilbert, 1995). In this problem, the magnetic field can be treated as passive. Furthermore, 
the viscosity-to-diffusivity ratio is often large enough for a sizable interval of scales between 
the viscous and the diffusive cut-offs to be present. That is the region of scales with the 
fastest growth rates of the magnetic fluctuations. Their dynamics, modeled by the passive 
advection of magnetic field by a large-scale (smooth) velocity field, will be described here. 

The dynamo process is caused by the stretching of fluid elements already extensively 
discussed above and the major new point to be noted is the role of the diffusion. In a perfect 
conductor, when the diffusion is absent, the magnetic field satisfies the same equation as the 
infinitesimal separation between two fluid particles (12): dB/dt = oB. Any chaotic flow 
would then produce dynamo, with the growth rate 



equal to the highest Lyapunov exponent Ai. Recall that the gradients of a scalar grow with 
the growth rate —A3 during the diffusionless stage. In fact, any real fluid has a nonzero 
diffusivity and, even though it can be very small, its effects may be dramatic. The long- 
standing problem solved by Chertkov, Falkovich et al. (1999) was whether the presence of 
a small, yet finite, diffusivity could stop the dynamo growth process at large times (as it is 
the case for the gradients of a scalar). 
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Our starting point is (99), expressing the magnetic field in terms of the stretching matrix 
W and the backward Lagrangian trajectory. In incompressible flow, matrices W and W are 
identically distributed and we do not distinguish them here. For example, the second-order 
correlation function is given by 



(%(r 12 ,t) = (B\ ri ,t) B*(r 2 ,t)) = (w il W^ m C l 2 m (R 12 (0;r 12 ,t),0)), (122) 



with the average taken over the velocity and the molecular noise. For the sake of simplicity, 
we assume that the initial statistics of B is homogeneous, Gaussian, of zero mean and of 
correlation length £. We concentrate on the behavior at scales r 12 <C £. For times less than 
td — l^l -1 l n (^/r d ), the Lagrangian separations R 12 <C £ and the magnetic field is stretched 
by the W matrix as in a perfect conductor, see (122). For longer times, the separation R 12 
can reach £, irrespective of its original value. This is the long-time asymptotic regime of 
interest, where the destinies of the scalar gradients and the magnetic field are different. 

It follows from (122) that the correlations are due to those realizations where Rn(0) < £. 
As in the Section III. A. 3, the initial separation r 12 should then be quasi-orthogonal to the 
expanding direction p 3 of W~ x and the fraction of solid angle occupied by those realizations is 
oc (£/ri 2 ) e p:i . Together with the e 2pi factor coming from the perfect conductor amplification, 
we thus obtain for the trace of the correlation function : 



with V(pi,p 2 ;t) as in (21). The integration is constrained by — \n(£/r 12 ) < p 2 , required 
for the separation along p 2 to remain smaller than £. Note that the gradients of a scalar 
field are stretched by the same W" 1 matrix that governs the growth of the Lagrangian 
separations. It is therefore impossible to increase the stretching factor of the gradient and 
keep the particle separation within the correlation length £ at the same time. That is 
why diffusion eventually kills all the gradients while the component B l that points into the 
direction of stretching survives and grows with V-B* perpendicular to it. This simple picture 
also explains the absence of dynamo in 2d incompressible flow, where the stretching in one 
direction necessarily means the contraction in the other one. 




(123) 
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Let us now consider the single-point moments (B 2n )(t). The 2n particles, all at the same 
point at time t, are split by the molecular diffusion by small separations of length 0(ra) 
in a time of the order td = r\l 'k near t. For the subsequent advection not to stretch the 
separations beyond £, the "diffusive" separations at time t — td should be quasi-orthogonal to 
the expanding direction p 3 . More exactly, they should form an angle < (^A"d) eP3 w hh the 
plane orthogonal to the expanding direction. Together with the pure conductor stretching 
factor, we are thus left with a contribution oc exp[n(2pi + p 3 )]. Two possible classes of 
Lagrangian trajectories should now be distinguished, depending on whether the angle formed 
by the "diffusive" separations with the p 2 direction is arbitrary or constrained to be small (see 
also Molchanov et al, 1985). For the former, the contribution is simply given by the average 
of the expression exp[n(pi — p 2 )] derived previously, with the constraint — \n£/r d < P2 
ensuring the control of the particle separation along p 2 . For the latter, the contribution is 
proportional to the average of exp[n(pi +2p 2 )]. Indeed, the condition of quasi-orthogonality 
to the p 2 direction contributes a np 2 term in the exponent and the remaining 2np 2 term 
is coming from the solenoidality condition V • B = 0. The magnetic field correlation is in 
fact proportional to the solenoidal projector and the component stretched by the W matrix, 
see (122), is C 2 X oc (1 — (V 2 ) _1 Vi)C(r). The realizations having the particle separations 
precisely aligned with the pi direction will therefore not contribute. For separations almost 
aligned to pi, one may show that the square of the angle with respect to p 1 appears in C^ 1 , 
thus giving the additional small factor (£/r d e P2 ) 2 . 

Which one of the two previous classes of Lagrangian trajectories dominates the moments 
depends on the specific form of the entropy function. For the growth rate (121), the situation 
is simpler as the average is dominated by the region around p« = A«. The average of the 
logarithm is indeed obtained by taking the limit n — > in (B 2n — l)/2n and the saddle 
point at large times sits at the minimum of the entropy function. The two previous classes 
of Lagrangian trajectories dominate for positive and negative A 2 , respectively. Using the 
identity — 0, we finally obtain 
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7 = min{(A 1 - A 2 ) /2, (A 2 - A 3 ) /2}. (124) 

The validity of this formula is restricted, first, by the condition that £exp(Ait) is still less 
than the viscous scale. The stretching in the third direction also imposes a restriction: 
the finiteness of the maximal possible size £ of the initial fluctuations gives the constraint 
£ exp(X 3 t) > r d . At larger time, (log|S| 2 ) decays. 

The most important conclusion coming from (124) is that the growth rate is always 
non-negative for a chaotic incompressible flow. Note that the growth rate vanishes if two 
of the Lyapunov exponents coincide, corresponding to the absence of dynamo for axially 
symmetric cases. For time-reversible and two-dimensional flows, the intermediate Lyapunov 
exponent vanishes and 7 = Ai/2. Note that 3d magnetic field does grow in a 2d flow; when, 
however, both the flow and the field are two-dimensional, one finds 7 = — Ai/2. For isotropic 
Navier-Stokes turbulence, numerical data suggest A 2 ~ Ai/4 (Girimaji and Pope, 1990) and 
the long-time growth rate is then 7 ps 3Ai/8. 

The moments of positive order all grow in a random incompressible flow with a nonzero 
Lyapunov exponent. Indeed, the curve E n = \n(B 2n )/2t is a convex function of n (due to 
Holder inequality) and it vanishes at the origin, where its derivative coincides by definition 
with the non-negative growth rate. Even when 7 = 0, the growth rates for n > are positive 
if the entropy function has a finite width. For n — 1 this was stated in (Gruzinov et ai, 
1996). As discussed previously, the behavior of the growth rate curve E n is nonuniversal 
and it depends on the specific form of the entropy function. For the Kraichnan case, we 
can use the result (32) for the entropy function and the calculation is elementary. The 
dominant contribution is coming from the average of exp[n(pi — p 2 )] and the p 2 integration 
is dominated by the lower bound — ln£/r^. The answer E 2 = 3Ai was first obtained by 
Kazantsev (1968). The general result is E n = \in(n + 4)/4, to be compared with the 
perfect conductor result Ain(2n + 3)/2. The difference between them formally means that 
the two limits of large times and small diffusivity do not commute (what is called "dissipative 
anomaly", see next Section). Multipoint correlation functions were calculated by Chertkov, 
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Falkovich et al. (1999). They reflect the prevailing strip structure of the magnetic field. 
An initially spherical blob evolves indeed into a strip structure, with the diffusive effects 
neutralizing one of the two directions that are contracted in a perfect conductor. The strips 
induce strong angular dependences and anomalous scalings similar to those described in 
Sect. III.B.3 below. 

7. Coil-stretch transition for polymer molecules in a random flow 

At equilibrium, a polymer molecule coils up into a spongy ball whose typical radius is 
kept at i?o by thermal noise. Being placed in a flow, such molecule is deformed into an 
elongated ellipsoid which can be characterized by its end-to-end extension R. As long as 
the elongation is much smaller than the total length of the molecule, the entropy is quadratic 
in R so that the molecule is brought back to its equilibrium shape by a damping linear in 
R. The equation for the elongation is as follows (Hinch, 1977) 

d t R+v VR= R Vv --R + rj. (125) 

T 

The left hand side describes advection of the molecule as a whole, the first term on the 
right hand side is responsible for stretching, r is the relaxation time and rj is the thermal 
noise with (rf (t) rf(0)) = S^S(t)R^/r. Since the size of the molecules is always much smaller 
than the viscous length then V« = a and one can solve (125) using the evolution matrix W 
introduced in Sect. II.B.l. At long enough times (when the initial condition is forgotten) the 
statistics of the elongation is given by R = / °° dsW(s)r](s) exp [—t/r\. We are interested 
in the tail of the PDF V(R) at R 3> Rq. The events contributing to it are related to the 
realizations with a long-time history of stretching where the variable p\ (corresponding to 
the largest Lyapunov exponent Ai) is large. The tail of the PDF is estimated analyzing the 
behavior of R/ Rq = / °° exp(pi(s) — s/r) ds. The realizations dominating the tail are those 
where pi(s) — s/r takes a sharp maximum at some time s* before relaxing to its typical 
negative values. The probability of those events is read from the large deviation expression 
(21): InV ~ — s*i?(s*~Vi( s *) ~~ where H is the entropy function. With logarithmic 
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accuracy one can then replace pi(s*) = \n(R/R ) + s*/r and what is left is just to find the 
maximum with respect to s*. The extremum value X* = (s*)" 1 \n(R/R ) is fixed by the 
saddle-point condition that H — X*H' should vanish at X* + r^ 1 — Ai. The final answer for 
the PDF is as follows 

V(R) oc R% R- l - a with a = H' (x* + r" 1 - A x ) . (126) 

The convexity of the entropy function ensures that a is positive if Ai < 1/r. 

In accordance with (126), the exponent a decreases when Ai increases and it tends to 
zero as Ai — > 1/r. In this region, the entropy function is quadratic and the exponent is 
expressed via the average value of p 1 and its dispersion only: a = 2(1 — Ait)/tA. The 
integral of the PDF diverges at large R as a tends to zero. The transition at Ai — > 1/r 
is called the coil-stretch transition as the majority of the polymer molecules got stretched. 
This stretching can be stopped by non-linear elastic effects or by the back reaction of the 
polymers onto the flow. The understanding of the coil-stretch transition goes back to the 
works by Lumley (1972, 1973). The power-law tail (126) has been derived by Balkovsky 
et al. (2000). The influence of non-linear effects on the statistics of the elongation was 
examined by Chertkov (2000). 

B. Cascades of a passive scalar 

This Section describes forced turbulence of a passive tracer statistically stationary in 
time and homogeneous in space. We consider the advection-diffusion equation 

d t 9 + (v ■ V)0 = nV 2 6 + (127) 

with the pumping ip assumed stationary, homogeneous, isotropic, Gaussian, of zero mean 
and with covariance 

(<p(r,t)<p(0,0)) = 6(t)$(r/e). (128) 
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The function $ is taken constant for r/£ < 1 and decaying rapidly for large ratios. The 
following considerations are valid for a pumping finite-correlated in time provided its corre- 
lation time in the Lagrangian frame is much smaller than the stretching time from a given 
scale to the pumping correlation scale I. Note that in most physical situations the sources do 
not move with the fluid so that the Lagrangian correlation time of the pumping is the min- 
imum between its Eulerian correlation time and £/V, where V is the typical fluid velocity. 
Most of the general features of the advection are however independent of the details of the 
pumping mechanism and its Gaussianity and ^-correlation are not a very serious restriction, 
as it will be shown in Section III.C.l. 

Equation (127) implies for incompressible velocities the balance relation for the "scalar 
energy" density e = 6> 2 /2: 

d t e + V- 3 = -e + 0, (129) 

where e = ^(V^) 2 is the rate of dissipation, = ip6 is that of energy injection, and j = 
\6 2 v — kO'VO is the flux density. In a steady state, the injection must be balanced by the 
diffusive dissipation, while the stretching and the contraction by the velocity provide for 
a steady cascade of the scalar from the pumping scale £ to the diffusion scale r d (where 
diffusion is comparable to advection). 

The advection-diffusion dynamics induces the Hopf equations of evolution for the equal- 
time correlation functions. For a white-in-time pumping, one obtains 

N 

dtfr ...9 N ) + ...v n -VJ n ...9 N ) 

71=1 

= k £ (e 1 . . . v% . . . e N ) + £ (e, . A . e N ) $„ m (130) 

n=l n,m n m 

in the shorthand notation 9 n = 9(r n ,t), $ nm = &(r nm /£) etc. These equations are clearly 
not closed since the left hand side involves the mixed correlators of the advected fields and 
the velocities. An exception is provided by the case of the Kraichnan ensemble of velocities 
where the mixed correlators may be expressed in terms of those containing only the advected 
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fields, see Sect. III.C.l below. The stationary version of the 2-point Hopf equation may be 
written in the form 



The relative strength of the two terms on the left hand side depends on the distance. For 
velocities scaling as A r v oc r a , the ratio of advection and diffusion terms Pe(r) = A r vr/n 
may be estimated as r a+l / k. In particular, Pe = Pe(£) is called the Peclet number, and 
the diffusion scale r d is defined by the relation Pe(r d ) = 1 . 

In the "diffusive interval" r i2 <C r d , the diffusion term dominates in the left hand side of 
(131). Taking the limit of vanishing separations, we infer that the mean dissipation rate is 
equal to the mean injection rate e = (k(V6>) 2 ) = |$(0). This illustrates the aforementioned 
phenomenon of the "dissipative anomaly" : the limit k — > of the mean dissipation rate is 
non-zero despite the explicit k factor in its definition. The "connective interval" r d <C ryi <C 
t widens up at increasing Peclet number. There, one may drop the diffusive term in (131) 
and thus obtain 



The expression (4) may be derived from the general flux relation (132) by the additional as- 
sumption of isotropy. The relation (132) states that the mean flux of 9 2 stays constant within 
the convective interval and expresses analytically the downscale scalar cascade. For the ve- 
locity scaling A r v oc r a , dimensional arguments suggest that A r 9 oc r^ 1-0 ^ 2 (Obukhov, 
1949; Corrsin, 1951). This relation gives a proper qualitative understanding that the degrees 
of roughness of the scalar and the velocity are complementary, yet it suggests a wrong scaling 
for the scalar structure functions of order higher than the second one, see Sect. III.C.l. 

Let us now derive the exact Lagrangian expressions for the scalar correlation functions. 
The scalar field along the Lagrangian trajectories R(t) changes as 



((t>i-Vi + t> 2 -V 2 )0i0 2 ) + 2«(v 1 0i-V 2 2 ) = $ 



(131) 




(132) 



d_ 
~dt 



6{R(t),t) = <p{R(t),t). 



(133) 
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The N-th order scalar correlation ( J] 9(r n ,t) ) = C N (r,t) is then given by 

\ P =i ' 

C w (c,f) = ( d Sl ... £<p(R N (s«),s«) ds N ), (134) 

with the Lagrangian trajectories satisfying the final conditions R n (t) = r n . For the sake 
of simplicity we have written down the expression for the case where the scalar field was 
vanishing at the initial time. If some of the distances among the particles get below the 
diffusive scale, the molecular noises in the Lagrangian trajectories become relevant and the 
averaging of (134) over their statistics is needed. 

The average over the Gaussian pumping in (134) gives for the pair correlation function: 

C 2 (rW) = Q\(R 12 (s)/£) ds^. (135) 

The function $ essentially restricts the integration to the time interval where R\ 2 is 
smaller than the injection length £. If the Lagrangian trajectories separate, the pair 
correlation reaches at long times the stationary form given by the same formula with 
t — oo. Simply speaking, the stationary pair correlation function of a tracer is pro- 
portional to the average time that two particles spent in the past within the correlation 
scale of the pumping (Falkovich and Lebedev, 1994). Similarly, the pair structure func- 
tion S 2 (r) = {(9i — 9 2 ) 2 ) = 2 [C 2 (0) — C 2 (r)\ is proportional to the time it takes for two 
coinciding particles to separate to a distance r. This is proportional to r 1_a for a scale- 
invariant velocity statistics with A r v oc r a , see (38), so that S 2 (r) is in agreement with the 
Obukhov-Corrsin dimensional prediction. 

Higher-order equal-time correlation functions are expressed similarly by using the Wick 
rule to average over the Gaussian forcing: 

C 2n (r,t) = (J*<S>(R 12 ( Sl ))d Sl ... j\{R {2n _ 1)2n {s n ))ds n ) + ... , (136) 

where the remaining average is over the velocity and the molecular noise ensemble and the 
dots stand for the other possible pairings of the 2n points. The correlation functions may 
be obtained from the generating functional 
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exp 



0(r,t)x(r) dr 



exp 



rfs / / $(R 12 (s)) x(ri)x(r2) dr x dr 2 . (137) 



They probe the statistics of times spent by fluid particles at distances smaller than I. 
In nonsmooth flows, the correlation functions at small scales, <C £, are dominated by 
the single-point contributions, corresponding to initially coinciding particles. This is the 
signature of the explosive separation of the trajectories. To pick up a strong dependence on 
the positions r, one has to study the structure functions which are determined by the time 
differences between different initial configurations. Conversely, the correlation functions 
at scales larger than £ are strongly dependent on the positions, as it will be shown in 
Sect. III.B.3. 



1. Passive scalar in a spatially smooth velocity 

In the rest of Section III.B, all the scales are supposed much smaller than the viscous 
scale of turbulence so that we may assume the velocity field to be spatially smooth and use 
the Lagrangian description developed in Sects. II. B and II. D. In the Batchelor regime, the 
backward evolution of the Lagrangian separation vector is given by #12(0) = W^(£) _1 ri 2 
(if we ignore diffusion) and it is dominated by the stretching rate pd at long times. The 

equation (135) takes then the asymptotic form 

t 

C 2 (r,t) » JdsJ<5>(e-^r/e)v( Pl ,...,p d ;s)dp l ...dp d . (138) 


The behavior of the interparticle distance crucially depends on the sign of Xj. For A^ < 0, 
the backward-in-time evolution separates the particles and leads in the limit t — > 00 to a 
well-defined steady state with the correlation function 

(6(t, 0) 9(t, r)) « 1^1-^(0) Hl/r) , (139) 

for r < t. This corresponds to the direct cascade. Conversely, if A^ > the particles 
contract and the pair correlation function grows proportionally to t. Note that the growing 
part is independent of r. This means that, in a flow contracting backwards in time, tracer 
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fluctuations grow at larger and larger scales, which is a signature of the so-called inverse 
cascade of a passive tracer. 

If the velocity ensemble is time-reversible, as it is the case for the 5-correlated model 
(57), then Aj = — \d-i+i and Ai and A^ have opposite sign. They will thus both change 
sign at the same value of the degree of compressibility p = d/A, see (60). This is peculiar 
for a short-correlated case and does not hold for an arbitrary velocity statistics. There, 
the change from stretching to contraction in the forward Lagrangian dynamics does not 
necessarily correspond to the change in the direction of the cascade for the passive tracer, 
related to the backward in time Lagrangian dynamics. 

2. Direct cascade, small scales 

We consider here the case X d < (that includes smooth incompressible flows) when 
the particles do separate backward in time and a steady state exists. We first treat 
the convective interval of distances between the diffusion scale r d and the pumping scale 
£. Deep inside the convective interval where r <C I, the statistics of the passive scalar 
tends to become Gaussian. Indeed, the reducible part in the 2n-point correlation function 
^$( e -pd(si) r 12 /£) ■ ■ ■ $(e^ d( * n) r {2n _i )2n /T)^, see (136) and (138), dominates the irreducible 
one for n <C n cr ps (IA^t,,) -1 \n(£/r). The reason is that the logarithmic factors are smaller 
for the irreducible than for the reducible contribution (Chertkov et ai, 1995a). The critical 
order n cr is given by the ratio between the time for the particles to separate from a typical 
distance r to £ and the correlation time t s of the stretching rate fluctuations. Since £ 3> r, 
the statistics of the passive tracer is Gaussian up to orders n cr ^> 1. For the single-point 
statistics, the scale appearing in the expression of the critical order n cr should be taken as 
r = r d . The structure functions are dominated by the forced solution rather than the zero 
modes in the convective interval: S 2 n = ([#(0) — 6{r)] 2n ) oc \n n (r/r d ) for n <^ \n(r/r d ). A 
complete expression (the forced solution plus the zero modes) for the four-point correlation 
function in the Kraichnan model can be found in Balkovsky, Chertkov et ai (1995). 
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Let us show now that the tails of the tracer PDF decay exponentially (Shraiman and 
Siggia, 1994; Chertkov et al, 1995a; Bernard et al, 1998; Balkovsky and Fouxon, 1999). The 
physical reasons behind this are transparent and most likely they apply also for a nonsmooth 
velocity. First, large values of the scalar can be achieved only if during a long interval of 
time the pumping works uninterrupted by stretching events that eventually bring diffusion 
into play. We are interested in the tails of the distribution, i.e. in intervals much longer than 
the typical stretching time from r d to I. Those rare events can be then considered as the 
result of a Poisson random process and the probability that no stretching occurs during an 
interval of length t is cx exp(— ci). Second, the values achieved by the scalar in those long 
intervals are Gaussian with variance $(0)t. Note that this is also valid for a non-Gaussian 
and finite-correlated pumping, provided t is larger than its correlation time. By integrating 
over the length of the no-stretching time intervals with the pumping-produced distribution 
of the scalar we finally obtain: V{9) oc / dt exp(-ct - 2 / 2 $(O)t) oc exp(-0^/2c/$(O)). 
This is valid for t < L 2 / k that is for 9 < ^/c$(0) L 2 / k. Interval of exponential behavior thus 
increases with the Peclet number. For a smooth case, the calculations have been carried 
out in detail and the result agrees with the previous arguments. Experimental data in 
(Jullien et al, 2000) confirm both the logarithmic form of the correlation functions and 
the exponential tails of the scalar PDF. In some experimental set-ups the aforementioned 
conditions for the exponential tails are not satisfied and a different behavior is observed, 
see for example (Jayesh and Warhaft, 1991). The physical reason is simple to grasp. The 
injection correlation time in those experiments is given by L/V, where L is the velocity 
integral scale and V is the typical velocity. The no-stretching times involved in the tail of 
the scalar distribution are of the order of W/V, where W is the width of the channel where 
the experiment is performed. Our previous arguments clearly require W ^> L. As the width 
of the channel is increased, the tails indeed tend to become exponential. 

For a 5-correlated strain, the calculation of the generating functional of the scalar cor- 
relators may be reduced to a quantum mechanical problem. In the Batchelor regime and 
in the limit of vanishing k , the exponent in the generating functional (137) may indeed be 
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rewritten as 

-J*V x (W(s)) ds = -\ J* ds J J <5>(W~\s)r 12 ) x (n) *(r 2 ) dr x dr 2 . (140) 

Recall that the matrices W form a stochastic process describing a diffusion on GL(d) 
(or on SL(d) in the incompressible case) with a generator M. The above formula may 
thus be interpreted as the Feynman-Kac expression for the integral J e t( -^~ Vx \l, W) dW 
of the heat kernel of M perturbed by the positive potential V x . As long as the trajectories 
separate backward in time, i.e. for < 0, the generating functional has a stationary limit, 
given by (140) with the time integral extending to infinity. The Feynman-Kac formula 
may be used to find the exponential rate of decay of the PDF V(9) oc e~ h \ e \ . As shown 
by Bernard et al. (1998), this involves the quantum-mechanical Hamiltonian — M — a 2 V x , 
where the positive operator — M has its spectrum starting at a strictly positive value and 
the negative potential tends to produce a bound state as the parameter a is increased. In 
the incompressible case, the decay rate b is characterized by the property that the ground 
state of the Hamiltonian has zero energy. For isotropic situations, the potential is only a 
function of the stretching rates of W and the quantum-mechanical problem reduces to the 
perturbation of the Calogero- Sutherland Hamiltonian by a potential, see Sect. II. B. 2. 

3. Direct cascade, large scales 

We consider here the scales r ^> £ in the steady state established under the condition 
\d < (Balkovsky et al, 1999b). From a general physical viewpoint, it is of interest to 
understand the properties of turbulence at scales larger than the pumping scale. A natural 
expectation is to have there an equilibrium equipartition with the effective temperature 
determined by the small-scale turbulence (Forster et al, 1977; Balkovsky, Falkovich et al, 
1995). The peculiarity of our problem is that we consider scalar fluctuations at scales larger 
than that of pumping yet smaller than the correlation length of the velocity field. This 
provides for an efficient mixing of the scalar even at those large scales. Although one can 
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find the simultaneous correlation functions of different orders, it is yet unclear if such a 
statistics can be described by any thermodynamical variational principle. 

The correlation functions of the scalar are proportional to the time spent by the La- 
grangian particles within the pumping scale. It follows that the statistics at r t is related 
to the probabilities of initially distant particles to come close. For spatially smooth random 
flow, such statistics turns out to be strongly intermittent and non-Gaussian. Another unex- 
pected feature in this regime is a total breakdown of scale invariance: not only the scaling 
exponents are anomalous and do not grow linearly with the order of the correlation function, 
but even fixed correlation functions are generally not scale invariant. The scaling exponents 
depend indeed on the angles between the vectors connecting the points. Note that the large- 
scale statistics of a scalar is scale- invariant in a nonsmooth velocity, see (Balkovsky et ai, 
1999b) and Sect.IILC.lv. 

What is the probability for the vector Ri 2 (t), that was once within the pumping corre- 
lation length £, to come exactly to the prescribed value r at time t? The advection makes 
a sphere of "pumping" volume £ d evolve into an elongated ellipsoid of the same volume. 
Ergodicity may be assumed provided that the stretching time A -1 ha(r/£) is larger than the 
strain correlation time. It follows that the probability for two points separated by r to be- 
long to a "piece" of scalar originated from the same pumping sphere behaves as the volume 
fraction (i/r) d . That gives the law of decrease of the two-point function: C 2 oc r~ d . 

The advection by spatially smooth velocities preserves straight lines. To determine the 
correlation functions of an arbitrary order when all the points lie on a line, it is enough to 
notice that the history of the stretching is the same for all the particles. Looking backward 
in time we may say that when the largest distance among the points was smaller than £, 
then all the other distances were as well. It follows that the correlation functions for a 
collinear geometry depend on the largest distance r among the points so that C 2n oc r~ d . 
This is true also when different pairs of points lie on parallel lines. Note that the exponent 
is n-independent which corresponds to a strong intermittency and an extreme anomalous 
scaling. The fact that C 2n 3> C 2 is due to the strong correlation of the points along their 
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common line. 

The opposite takes place for non-collinear geometries, namely the stretching of different 
non-parallel vectors is generally anti-correlated in the incompressible case due to the volume 
conservation. The ci-volume e^...^ RY 2 . . . F£y d is indeed preserved for (d + 1) Lagrangian 
trajectories R n {t) and, for d = 2 and any three trajectories, the area R\ 2 R{ 3 of the 
triangle defined by the three particles remains constant. The anti-correlation due to the area 
conservation may then be easily understood and the scaling for non-collinear geometries at 
d = 2 may be determined. Since the area of any triangle is conserved, three points that 
form a triangle with area A ^> £ 2 will never come within the pumping correlation length. 
In the presence of a triple correlator $ 3 for a non Gaussian 5-correlated pumping, the triple 
correlation function of a scalar 



is determined by the asymptotic behavior of <3> 3 at 3> I. For example, if $ 3 has a 
Gaussian tail, then C 3 oc exp(— A/£ 2 ). On the other hand, the correlation functions decrease 
as r~ 2 for a collinear geometry. We conclude that C 3 as a function of the angle between the 
vectors r 12 and r 13 has a sharp maximum at zero and rapidly decreases within an interval 
of width of the order £ 2 /r 2 <C 1. 

Similar considerations apply for the fourth-order correlation function. Note that, unlike 
for the 3-point function, there are now reducible contributions. Consider, for instance, that 



coming from ( / / $(i? 12 (si)) $(-R 3 4(s 2 )) dsids 2 ). Since the area of the polygon defined by 
x o o / 

the four particles is conserved throughout the evolution, the answer is again crucially de- 
pendent on the relation between the area and £ 2 . The events contributing to the correlation 
function C 4 are those where, during the evolution, i? 12 became of the order £ and then, at 
some other moment of time, _R 34 reached I. The probability for the first event to happen is 
£ 2 jr\ 2 . When this happens, the area preservation makes i? 34 ~ ri 2 r 34 /£. The probability 
for this separation to subsequently reduce to £ is oc £ 4 /rf 2 r u- The total probability can be 
thus estimated as £ 6 /r 6 , where r is the typical value of the separations r\y Remark that the 




(141) 
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naive Gaussian estimation £ 4 /r 4 is much smaller than the collinear answer and yet much 
larger than the non-collinear one. 

The previous arguments can be readily extended to an arbitrary number of non-collinear 
pairs. In accordance with (136), the realizations contributing to the correlation are those 
where the separations reduce down to t during the evolution process. Suppose that this 
happens first for Ri 2 . Such process was already explained in the consideration of the pair 
correlation function and occurs with probability (£/ri 2 ) 2 . All the remaining separations will 
then be larger than their initial values by a factor r\ijl, due to the conservation law of the 
triangular areas. Next, we should reduce, say, -R34 from r 34 ri 2 /£ down to the integral scale 
i. Such a process occurs with probability (^ 2 /(ri 2 r 34 )) 2 . When this happens, all the other 
separations are larger than their initial values by a factor r u /£. Repeating the process, 
we come to the final answer C 2 „ oc (£/r) 4 ™ -2 , where r is again the typical value of the 
separations r\y 

The above analysis is easy to generalize for arbitrary geometries. The points are divided 
into sets consisting of the pairs of points with parallel separations r ^ (more precisely, forming 
angles smaller than £ 2 /r 2 ). The points within a given set behave as a single separation during 
the Lagrangian evolution. The order n in the previous formulae should then be replaced by 
the (minimal) number of sets. The estimates obtained above are supported by the rigorous 
calculations in (Balkovsky et al, 1999b). 

In 2d, the area conservation allowed to get the scaling without calculations. This is 
related to the fact that there is a single Lyapunov exponent. When d > 2, we have only 
the conservation of ci-dimensional volumes and hence more freedom in the dynamics. For 
example, the area of a triangle can change during the evolution and the three-point cor- 
relation function for a non-collinear geometry is not necessarily suppressed. Nevertheless, 
the anti-correlation between different Lagrangian trajectories is still present and the 3-point 
exponent is expected to be larger than the naive estimate 2d. The answer for the Kraich- 
nan model d + (d — l)\Jd/ {d — 2) is determined by the whole hierarchy of the Lyapunov 
exponents (Balkovsky et al, 1999b). In the limit of large dimensions, the anti-correlation 
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tends to disappear and the answer approaches 2d. The four-point correlation function is 
also determined by the joint evolution of two distances, which results for large d in the same 
value of the exponent. 

To conclude this section, we briefly comment on the case A^ > when the particles 
approach rather than separate backward in time. Here, an inversion of what has been 
described for the direct cascade takes place: the scalar correlation functions are logarithmic 
and the PDF has a wide Gaussian core at r > £, while the statistics is strongly non-Gaussian 
at small scales (Chertkov et al, 1998). Since the scalar fluctuations injected at £ propagate 
upscale, small-scale diffusion is negligible and some large-scale damping (say, by friction) is 
needed to provide for a steady state, see also Sect. III.E below. 

4. Statistics of the dissipation 

We now describe the PDF's of the scalar gradients u> = V# and of the dissipation e = k uo 2 
in the steady state of a direct cascade arising under the action of a large-scale pumping. 
We consider a smooth velocity field, i.e. both the Schmidt/Prandtl and the Peclet numbers 
are assumed large. As we remarked in Sect. III. A. 5, the scalar gradients can be estimated 
as 9e~ Pd /£, where 9 is the scalar value and e Pd £ is the smallest (diffusive) scale. The tails 
of the gradient PDF are controlled by the large values of 9 and — pd- The statistics of the 
former depends both on the pumping and the velocity and that of the latter only on the 
velocity. The key remark for solving the problem was made by (Chertkov et al, 1997 and 
1998): since 9 and — p d fluctuate on very separated time scales ( |A rf | _1 \n(£/r d ) and A^ 1 , 
respectively), their fluctuations may be analyzed separately. The PDF of the scalar has been 
shown in Sect. III.B.2 to decay exponentially. On the other hand, large negative values of 
Pd are determined by the tail of its stationary distribution, see (27). For a Gaussian short- 
correlated strain this tail is oc exp[— const. e~ 2pd ] (Chertkov et al, 1998) and the moments 
of the gradients are (cu n ) oc (9 n ) {exp(-npd)) oc n 3ri//2 . The ensuing behavior (e n ) oc n 3n 
corresponds to a stretched-exponential tail for the PDF of the dissipation 
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\nV(e) oc -e 1/3 . (142) 

The detailed calculation for the Kraichnan model as well as a comparison with numerical and 
experimental data can be found in (Chertkov, Kolokolov et al, 1998; Chertkov, Falkovich 
et al, 1998; Gamba & Kolokolov, 1998). The general case of a smooth flow with arbitrary 
statistics was considered in (Balkovsky and Fouxon, 1999). 

It is instructive to compare the stretched-exponential PDF of the gradients in a steady 
state with the lognormal PDF described in Sect. III. A. 5 for the initial diffusionless growth. 
Intermittency builds up during the initial stage, i.e. the higher the moment, the faster it 
grows. On the other hand, the higher the moment, the shorter is the breakdown time of the 
diffusionless approximation. This time behaves for example as (n + 2) _1 in the Kraichnan 
model. Since higher moments stop growing earlier than lower ones, the tails of the PDF 
become steeper and the intermittency is weaker in the steady state. 

C. Passive fields in the inertial interval of turbulence 

For smooth velocities, the single-point statistics of the advected quantities could be 
inferred from the knowledge of the stretching rates characterizing the Lagrangian flow in 
the infinitesimal neighborhood of a fixed trajectory. This was also true for the multipoint 
statistics as long as all the scales involved were smaller than the viscous scale of the velocity, 
i.e. in the Batchelor regime. In this Section we shall analyze advection phenomena, mostly 
of scalars, in the inertial interval of scales where the velocities become effectively nonsmooth. 
As discussed in Sect. II. C, the explosive separation of the trajectories in nonsmooth velocities 
blows up interparticle separations from infinitesimal to finite values in a finite time. This 
phenomenon plays an essential role in maintaining the dissipation of conserved quantities 
nonzero even when the diffusivity k — > . The statistics of the advected fields is consequently 
more difficult to analyze, as we discuss below. 
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1. Passive scalar in the Kraichnan model 

The Kraichnan ensemble of Gaussian white-in-time velocities permits an exact analy- 
sis of the nonsmooth case and a deeper insight into subtle features of the advection, like 
intermittency and anomalous scaling. Those aspects are directly related to the collective 
behavior of the particle trajectories studied in the first part of the review. Important lessons 
learned from the model will be discussed in next Sections in a more general context, 
i). Hopf equations. The simplifying feature associated to the Kraichnan velocities is a 
reduction of the corresponding Hopf equations to a closed recursive system involving only 
correlators of the advected fields. This is due to the temporal decorrelation of the velocity 
and the ensuing Markov property of the Lagrangian trajectories. Let us consider for example 
the evolution equation (127) for a scalar field. For the Kraichnan model, it becomes a 
stochastic differential equation. As mentioned in Sect. II. B. 2, one may view white-in-time 
velocities as the scaling limit of ensembles with short time correlations. The very fact 
that v(t)dt tends to become of the order (dt) 1 ^ 2 calls for a regularization. For velocity 
ensembles invariant under time reversal, the relevant convention is that of Stratonovich (see 
the Appendix). Interpreting (127) within this convention and applying the rules of stochastic 
differential calculus, one obtains the equation for the scalar correlation functions: 

d t C N {r) = M N C N (r) + ^C^K^^rJ^Mf). (143) 

n<m n m 

Here, the differential operator A4 N is the same 9 as in (69) and it may be formally obtained 
from the second term on the left hand side of (130) by a Gaussian integration by parts. Note 
the absence of any closure problem for the triangular system of equations (143): once the 
lower-point functions have been found, the iV-point correlation function satisfies a closed 
equation. For spatially homogeneous situations, the operators A4 N may be replaced by 
their restrictions Ai M to the translation-invariant sector. It follows from their definition 



9 The n = m terms would drop out of the expression for M N in the ltd convention for (127). 
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(70) and the velocity correlation function (48) that the equations (143) are then invariant 
with respect to the rescalings 

r^Xr, £^X£, t — > A 2_? t , k -> A f k , 9 ^ X~^ 9 . (144) 

This straightforward observation implies scaling relations between the stationary correlators: 
C N (Xr; X^k, X£) = At C N (r; k,£). 

ii). Pair correlator. For the isotropic pair correlation function, the equation (143) takes 
the form: 

d t C 2 {r) -r l ~ d d r [(d- l)D ir d - 1+ t + 2Kr d ~ 1 ] d r C 2 (r) = <$>{r/£) , (145) 

see (53). The ratio of the advective and the diffusive terms is of order unity at the 
diffusion scale = [2k/ (d — 1)D\\ 1 ^. For the Kraichnan model, the Peclet number 
Pe = (d— 1)Di£^/2k > 1 as we assume the scale of pumping much larger than that of 
diffusion. The stationary form of (145) becomes an ordinary differential equation (Kraich- 
nan, 1968) that may be easily integrated with the two boundary conditions of zero at infinity 
and finiteness at the origin: 

c *w = <5=k / irdr / y"- 1 "y ■ < 146 > 

r X ^ ' d o 

Even without knowledge of the explicit form (146), it is easy to draw from (145) general 
conclusions, as for time-correlated velocities. Taking the limit r — > for k > 0,we infer 
the mean scalar energy balance 

<9*e + e = $(0)/2, (147) 

where e = (8 2 ) /2 and e = (/t(V#) 2 ). In the stationary state, the dissipation balances the 
injection. On the other hand, for r ^> (or for any r > in the k — > limit) we may drop 
the diffusive term in the Hopf equation (145). For r <C £, we thus obtain the Kraichnan 
model formulation of the Yaglom relation (132) expressing the constancy of the downscale 
flux: 
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-{d-l)D 1 -^- I d r r d - 1+ ^d r C 2 (r,t) « $(0). (148) 

To obtain the balance relation (147) for vanishing k as the r — > limit of (145), one has 
to define the limiting dissipation field by the operator product expansion 

lim K{yef(r) = \ lim d ij (r - r') V^(r) V Mr') . (149) 

The relation (149), encoding the dissipative anomaly, holds in general correlation functions 
away from other insertions (Bernard et al, 1996). 

Let us discuss now the solution (146) in more detail. There are three intervals of distinct 
behavior. First, at large scales r £, the pair correlation function is given by 

where $ = / °° dy. This may be thought of as the Rayleigh- Jeans equipartition 

(9(k) 6{k')) = S(k + k') ®£ d /£l(k) with Q(k) oc k 2 ~^ and the temperature proportional to 
<t£ d . Note that the right hand side of (150) is a zero mode of M 2 away from the origin: 
r 1 ~ d d r r d ~ 1+ ^d r r 2 ~^ d oc S(r). Second, in the convective interval r& <C r <C £, the pair 
correlator is equal to a constant (the genuine zero mode of A4 2 ) plus an inhomogeneous 
part: 

C 2 (r) « A 2 £^ - (2 _^_ 1)Di mr 2 -*, (151) 

where A 2 = $(0)/(2 — + £ — 2)(d — l)Di. The leading constant term drops out of the 
structure function: 

S 2 (r) = 2 [C 2 (0) - 2C 2 (r)} » (2 _^ ( l 1)Dl $(0) r 2 -« . (152) 

Note that the last expression is independent of both k and I and it depends on the pumping 
through the mean injection rate only, i.e. S 2 (r) is universal. Its scaling exponent ( 2 = 2 — 1; 
is fixed by the dimensional rescaling properties (144) or, equivalently, by the scaling of the 
separation time of the Lagrangian trajectories, see (55). As remarked before, the degrees of 
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roughness of the scalar and the velocity turn out to be complementary: a smooth velocity 
corresponds to a rough scalar and vice versa. Finally, in the diffusive interval r <C r d , 
the pair correlation function is dominated by a constant and the structure function S^r) ~ 
_L_ $(Q) r 2 . Note that S 2 (r) is not analytic at the origin, though. Its expansion in r contains 
noninteger powers of order higher than the second, due to the non-smoothness of the velocity 
down to the smallest scales. The analyticity is recovered if we keep a finite viscous cutoff 
for the velocity. 

In the limit k — > 0, the diffusive interval disappears and the pair correlator is given by 
(146) with the diffusion scale set to zero. The mean square of the scalar 6*2(0) remains 
finite for finite £ but diverges in the £ — > 00 limit that exists only for the structure function. 
Recall from Sect. III. B.l that C 2 (r) has the interpretation of the mean time that two La- 
grangian trajectories take to separate from distance r to t. The finite value of the correlation 
function at the origin is therefore another manifestation of the explosive separation of the 
Lagrangian trajectories. 

The solution (146) for the pair correlator and most of the above discussion remain valid 
also for £ = 2, i.e. for smooth Kraichnan velocities. A notable difference should be stressed, 
though. In smooth velocities and for k — > 0, the mean time of separation of two Lagrangian 
trajectories diverges logarithmically as their initial distance tends to vanish. The pair corre- 
lator has a logarithmic divergence at the origin, implying that (9 2 ) is infinite in the stationary 
state with k = 0. Indeed, it is the <9iC 2 (0) term that balances the right hand side of (145) at 
finite times and r = 0. As the diffusivity vanishes, the variance (8 2 ) keeps growing linearly 
in time with the rate $(0) and the mean dissipation tends to zero: no dissipative anomaly 
is present at finite times. The anomaly occurs only in the stationary state that takes longer 
and longer to achieve for smaller r. In mathematical terms: lim lim e 7^ lim lim e, with 
the left hand side vanishing and the right hand side equal to the mean injection rate. The 
physics behind this difference is clear. The dissipative anomaly for nonsmooth velocities is 
due to the non-uniqueness of the Lagrangian trajectories, see Sect. II. C. The incompressible 
version of (100) implies that, in the absence of forcing and diffusion, 
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J 6 2 (r', 0) dr' - J 6 2 (r', t) dr' = j dr J p(r, t;R,0\ v) [e(R, 0) - 0(r, t)] 2 dR > . (153) 

The equality holds if and only if, for almost all r, the scalar is constant on the support of the 
measure p(r, t; R, | v) dR giving the distribution of the initial positions of the Lagrangian 
trajectories ending at r at time t. In plain language, / 9 2 dr is conserved if and only if the 
Lagrangian trajectories are uniquely determined by the final condition. This is the case for 
smooth velocities and no dissipation takes place for k = as long as / 6 2 dr is finite. When 
the latter becomes infinite (as in the stationary state), the above inequalities become void 
and the dissipation may persist in the limit of vanishing k even for a smooth flow. 

For £ = 0, the equation (146) still gives the stationary pair correlation function if d > 3. 
The distinction between the behavior in the convective and the diffusive regimes disappears. 
The overall behavior becomes diffusive with the stationary equal-time correlation functions 
coinciding with those of the forced diffusion d t 9 = (\(d — \)Di + /t)V 2 # + tp (Gaw§dzki 
and Kupiainen, 1996). In d — 2, the pair correlation function has a constant contribu- 
tion growing logarithmically in time but the structure function does stabilize, as in forced 
diffusion. 

iii). Higher correlators and zero modes. Let us consider the evolution of higher-order 
scalar correlation functions C N , assumed to decay rapidly in the space variables at the initial 
time. At long times, the correlation functions will then approach a stationary form given by 
the recursive relation 

C N (r) = I G N (r,R) £ C N _ 2 {R,, R N ) <$>{R n JZ) dR (154) 



for even N and vanishing for odd N (by the 9 \— > —9 symmetry). Here, G N = j e n dt 

o 

are the operators inverse to — M. N . The above formulae give specific solutions of the station- 
ary Hopf equations (143) that, alone, determine solutions only up to zero modes of operators 
M N . We are interested in the scaling properties of the stationary correlation function C N 
in the convective interval. If the correlation functions were becoming independent of k and 
i in this interval (mathematically, if the limits k — > and £ — > oo of the functions existed), 
the scaling behavior would follow from the dimensional relation (144): 
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C N (Xr) = \-^C N {r). (155) 

This would be the Kraichnan-model version of the normal Kolmogorov-Obukhov-Corrsin 
scaling (Obukhov, 1949; Corrsin, 1951). The k — > limit of the stationary correlation 
functions does exist and the ^-dependence drops out of the expressions in the convective 
interval, as for the pair correlator. The limit is given by the formulae (154) with the k = 
versions of G N (Hakulinen, 2000). Note in passing that the advection preserves any power 
of the scalar so that dissipative anomalies are present also for orders higher than the second. 
The existence of the zero diffusivity limit means that possible violations of the normal 
scaling in the convective interval may only come from a singularity of the limit £ — > oo. 
In fact, this was already the case for C2, dominated by the constant term that diverged 
as £ increases. The constant dropped out, however, from the pair structure function (152) 
that did not depend on £ and, consequently, scaled dimensionally. Concerning higher-order 
scalar structure functions, Kraichnan (1994) was the first to argue in favor of their anomalous 
scaling. His paper steered a renewed interest in the problem which led to the discovery by 
Chertkov et al. (1995b), Gaw§dzki and Kupiainen (1995) and Shraiman and Siggia (1995) 
of a simple mechanism to avoid normal scaling: the domination of the correlation functions 
by scaling zero modes of the operators A4 N . For small £ , Gaw§dzki and Kupiainen (1995) 
and Bernard et al. (1996) showed that in the convective interval 

C N (r) = A N £ A »f Nfi (r) + C' N {r) + o{£) + [...]. (156) 

Above, f N0 is the irreducible isotropic zero mode of scaling dimension ( N0 = ^(2 — £) — A^, 
see Sect II. E. 3, the term C' N is scaling with the normal dimension y(2 — £), and [. . .] 
stands for reducible contributions depending only on a subset of points. The anomalous 
corrections A^ = N 2 ^^ £ + 0{^ 2 ) are positive for small £ , see (84). A similar result A^ = 
jV( ^~ 2 - ) + O(Js) was established by Chertkov et al. (1995b) and Chertkov and Falkovich 
(1996) for large space dimensionalities. For A^ > 0, the first term on the right hand side of 
(156) is dominating the second one for large £ or, equivalently, at short distances for k — 0. 
The analytic origin of the zero- mode dominance of the stationary correlation functions (154) 
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lies in the asymptotic short-distance expansion (79) of the kernels of G N (Bernard et ai, 
1998). The dominant zero mode f N0 is the irreducible term in (79) with the lowest scaling 
dimension. The reducible terms [ . . . ] drop out of the correlators of scalar differences, e.g. 
in the iV-point structure functions. The latter are dominated by the contribution from f N : 

S N (r) = ([6{r) - 6{0)] N ) oc £ A " r c " , (157) 

with Cjv = Cjvo- The physical meaning of zero- mode dominance is transparent. Any structure 
function is a difference between the terms with different number of particles coming at the 
points 1 and 2 at time t, like, for instance, S 3 = 3(6l6 2 — 6i6%). Under the (backward-in-time) 
Lagrangian evolution, this difference decreases as (r/£)^ N -° because of shape relaxation with 
the slowest term due to the irreducible zero mode. The structure function is thus given by 
the total temporal factor £ N ( 2 -Z)/ 2 multiplied by (r/£) CN >°. 

The scaling exponents in (157) are universal in the sense that they do not depend on 
the shape of the pumping correlation functions $(r). The coefficients A N in (156), as well 
as the proportionality constants in (157), are, however, non-universal. Numerical analysis, 
see Sect. III.D.2, indicates that the previous framework applies for all < £ < 2 at any 
space dimensionality, with the anomalous corrections continuing to be strictly positive 
for N > 2. That implies the small-scale intermittency of the scalar field: the ratios S 2n / S% 
grow as r decreases. At orders N 3> (2 — the scaling exponents (n tend to saturate 
to a constant, see Sects. III.C.2 and III.D.2 below. 

In many practical situations the scalar is forced in an anisotropic way. Shraiman and 
Siggia (1994, 1995) have proposed a simple way to account for the anisotropy. They sub- 
tracted from the scalar field an anisotropic background by defining 9'{r) = 9{r) — g ■ r , 
with g a fixed vector. It follows from the unforced equation (92) that 

d t 6' + v -V6' - kV 2 6' = -g -v, (158) 

with the term on the right hand side giving the effective pumping. In Kraichnan velocities, 
the translation invariance of the equal-time correlators of 9' is preserved by the evolution 
with the Hopf equations taking the form 
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d t C N (r) = M N C N (r) + 2 ]T C^fo, . x r w ) g i9j D^(r nm ) 

-Eft ^W^^K,.,,^) (159) 

ra,m n 

in the homogeneous sector. The stationary correlation functions of 9' which arise at long 
times if the initial correlation functions decay in the space variables, may be analyzed as 
before. In the absence of the 9' \— > —9' symmetry, the odd correlators are not anymore 
constrained to be zero. Still, the stationary 1-point function vanishes so that the scalar 
mean is preimposed: (9(r)) = g r. For the 2-point function, the solution remains the same 
as in the isotropic case, with the forcing correlation function simply replaced by 2 g^j D^{r) 
and approximately equal to the constant 2D g 2 in the convective interval. The 3-point 
function is 

C 3 (r) = - / G 3 (r,R) ]T 9i d^(R nm ) V^C 2 ( J R 1 , R 3 ) dR. (160) 

n,m n 

The dimensional scaling would imply that C 3 (Ar) = A 3 ~ ? C 3 (r) in the convective interval 
since VC 2 (r) scales there as r 1 ^. Instead, for £ close to 2, the 3-point function is 
dominated by the angular momentum j = 1 zero mode of with scaling dimension 
2 + o(2 — £), as shown by Pumir et al. (1997). A similar picture arises from the perturbative 
analysis around £ = (Pumir, 1996 and 1997), around d = oo (Gutman and Balkovsky, 
1996), and from the numerical study of the whole interval of £ values for d = 2 and 
d = 3 (Pumir, 1997), see Sect. II. E. 5. As will be discussed in Section III.F, the zero mode 
mechanism is likely to be responsible for the experimentally observed persistence of the 
anisotropies, see e.g. (Warhaft, 2000). 

It is instructive to analyze the limiting cases £ = 0,2 and d = oo from the viewpoint 
of the statistics of the scalar. Since the field at any point is the superposition of contribu- 
tions brought from d directions, it follows from the Central Limit Theorem that the scalar 
statistics becomes Gaussian as the space dimensionality d increases. In the case £ = 0, an 
irregular velocity field acts like Brownian motion. The corresponding turbulent transport 
process is normal diffusion and the Gaussianity of the scalar statistics follows from that of 
the input. What is general in the both previous limits is that the degree of Gaussianity, 
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as measured, say, by the flatness S4/S2, is scale-independent. Conversely, we have seen in 
Sect. III.B.l that ln(£/r) is the parameter of Gaussianity in the Batchelor limit with the 
statistics becoming Gaussian at small scales whatever the input statistics. The key here is in 
the temporal rather than the spatial behavior. Since the stretching in a smooth velocity field 
is exponential, the cascade time is growing logarithmically as the scale decreases. That leads 
to the essential difference: at small yet nonzero £/d the degree of non- Gaussianity increases 
downscales, while at small (2 — £) it first decreases downscales until ln(£/r) xs 1/(2 — £) 
and then it starts to increase. Note that the interval of decrease grows as the Batchelor 
regime is approached. Already that simple reasoning suggests that the perturbation theory 
is singular in the limit £ — > 2, which is formally manifested in the quasi- singularities of the 
many-point correlation functions for collinear geometries (Balkovsky, Chertkov et al, 1995). 

The anomalous exponents determine also the moments of the dissipation field e = 
k(V6) 2 . A straightforward analysis of (154) indicates that (e n ) = c n e a (£/r d ) A2n (Chertkov 
et al, 1995b; Chertkov and Falkovich, 1996), where e is the mean dissipation rate. The 
dimensionless constants c n are determined by the fluctuations of the dissipation scale and, 
most likely, they are of the form n qn with yet unknown q. In the perturbative domain 
n <C (2 — £)(£/£, the anomalies A 2n are a quadratic function of the order and the corre- 
sponding part of the dissipation PDF is close to lognormal (Chertkov and Falkovich, 1996). 
The form of the distant tails of the PDF are still unknown. 

iv). Operator product expansion. While the irreducible zero mode dominates the 
respective structure function, all the zero modes may be naturally incorporated into an 
operator product expansion (OPE) of the scalar correlation functions. There has been many 
attempts to use this powerful tool of quantum field theory (Wilson, 1969) in the context of 
turbulence, see (Eyink, 1993; Adzhemyan et al, 1996; Polyakov, 1993 and 1995). We briefly 
describe here a general direction for accomplishing that for the problem of scalar advection 
(Chertkov and Falkovich, 1996; Adzhemyan et al, 1998; Zamolodchikov et al, 2000). Let 
{O a } be a set of local observables (which contains all spatial derivatives of any field already 
included) . The existence of OPE presumes that 
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O a (r)O b (r') = Y,C: b (r-r')O c (r>), (161) 

k 

which is understood as the following relations among the correlation functions 

(O a (r) O b (r>) ...) = £ C: b (r - r>) (O c (r>) . . .) . (162) 

k 

The sum represents the correlation function in the left-hand side if \r — r'\ is small enough. 
Renormalization symmetry if — > Aip, 9 — > A9 allows one to classify the operators (fluctuat- 
ing fields) by degrees: has degree n if — > A™ 0( n ) under the transformation. The 
OPE conserves the degree and is supposed to be scale invariant in the convective interval. 
This means that one may choose a basis of the observables in such a way that O a has "di- 
mension" d a , and the OPE is invariant under the transformation O a (r) — > \ da O a (\r) so 
that its coefficient functions scale: C^ b (X(r — r')) = \ dc ^ da ^ db C^ b (r — r') . Besides, functions 
C c ah are supposed to be pumping-independent with the whole dependence on pumping car- 
ried by the expectation values (O c ) oc l~ dc . The scale invariance can be ("spontaneously") 
broken at the level of correlation functions if some of the fields with nonzero dimension 
develop nonzero expectation values. The dimension of 6 N is AT(£ — 2)/2. The operators 
can be organized into strings, each with the primary operator O a with the lowest dimen- 
sion d a and its descendants with the dimensionalities d a + n{2 — £). A natural conjecture 
is that there is one-to-one correspondence between the primary operators of degree N and 
the zero modes of M N . The dimensions of such primaries are minus the anomalous di- 
mensions A's of the zero modes f a and are therefore negative. By fusing — 1 times one 
gets 6 1 . . .6 N = fa(r 1 , . . . , r N ) O a (r JV ) + . . . , where the dots include the derivatives and 
descendants of Q a . 

For n = 2, one has only one primary field 8 2 , and its descendents | V m 9\ 2 . For n = 4, there 
is an infinity of primaries. Only 8 A , e 2 and e6'V 2 9 — d 2 e/(d 2 + 2) have nonzero expectation 
values. The operators with zero expectation values correspond to the operators with more 
derivatives than twice the degree (that is with the order of the angular harmonic in the 
respective zero mode being larger than the number of particles, in terms of Sect. II. E). 
Building an OPE explicitly and identifying its algebraic nature remains a task for the future. 
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v) . Large scales. The scalar correlation functions at scales larger than that of the pumping 
decay by power laws. The pair correlation function is given by (150). Recall that applying 
M. 2 on it, we obtain a contact term oc S(r). Concerning higher-order correlation functions, 
straight lines are not preserved in a nonsmooth flow and no strong angular dependencies 
of the type encountered in the smooth case are thus expected. To determine the scaling 
behavior of the correlation functions, it is therefore enough to focus on a specific geometry. 
Consider for instance the equation M^C^r) = X x( r 'ij) C 2 (rki) for the fourth order func- 
tion. A convenient geometry to analyze is that with one distance among the points, say ri 2 , 
much smaller than the other r±j, whose typical value is R. At the dominant order in r^/R, 
the solution of the equation is C 4 oc C2(r 12 )C 2 (R) ~ (r 12 R) 2 ~^~ d . Similar arguments apply 
to arbitrary orders. We conclude that the scalar statistics at r 3> I is scale- invariant, i.e. 
C2n(Ar) = A"( 2 ^~ d ) C 2n (r) as A — > oo. Note that the statistics is generally non-Gaussian 
when the distances between the points are comparable. As £ increases from zero to two, the 
deviations from the Gaussianity starts from zero and reach their maximum for the smooth 
case described in Sect. III.B.3. 

vi) . Non-Gaussian and finite correlated pumping. The fact that the scalar corre- 
lation functions in the convective interval are dominated by zero modes indicates that the 
hypotheses of Gaussianity and ^-correlation of the pumping are not crucial. The purpose of 
this Section is to give some more details on how they might be relaxed. The new point to be 
taken into account is that the pumping has now a finite correlation time r p and irreducible 
contributions are present. The situation with the second order is quite simple. The injection 
rate of 9 2 is ((f) 9) and its value defines the mean dissipation rate e at the stationary state. 
The only difference is that its value cannot be estimated a priori as $(0). Let us then 
consider the behavior at higher orders, whose typical example is the fourth. Its general flux 
relation, derived similarly as (131), reads 



(»i • Vi + v 2 ■ v 2 ) elel) + K (0102 [vj + vj|0i0 2 ) = MA + v20 2 o\) . (163) 

Taking the limit of coinciding points, we get the production rate of 4 . It involves the usual 
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reducible contribution 3eC2(0) and an irreducible one. The ratio of the two is estimated 
as C 2 (0)/r p . The non-Gaussianity of the pumping is irrelevant as long as r p is smaller 
than the time for the particles to separate from the diffusive to the integral scale. The 
smooth and nonsmooth cases need to be distinguished. For the former, the separation 
time is logarithmically large and the previous condition is always satisfied. Indeed, the 
reducible part of the injection rate in (163) necessarily contains 2(6*16*2) [(^1^2) + (^26*1)] — 
2e 2 A^~ 1 ln(£/ri 2 ). Since the correlation function grows as r 12 decreases, one can always neglect 
the constant irreducible contribution for small enough separations. Similarly, the input rate 
of all even moments up to iV ~ ln(£/r^) is determined by e. The fact that the fluxes of higher 
integrals are not constant in the convective interval was called the effect of "distributed 
pumping" in (Falkovich, 1994; Falkovich and Lebedev, 1994). 

In the nonsmooth case, the cascade time is finite and the irreducible contributions might 
be relevant. They affect the statistics of the scalar yet, of course, not the scaling of the zero 
modes. The fourth order correlation function C 4 acquires for example extra terms propor- 
tional to J2 r ij ■ They contribute to the fourth order cumulant but not to the structure 
function S4. The existence of those extra terms in the correlation function affects the match- 
ing conditions at the pumping scale, though. We conclude that the numerical coefficients A N 
in the structure functions S N = A N r^ £ a n generally depend on all the irreducible pumping 
contributions of order m < N . 

2. Instanton formalism for the Kraichnan model 

Since the perturbative approaches in Section II. E. 5 are all limited to finite orders, it 
is natural to look for alternative methods to capture the scaling exponents in the non- 
perturbative domain N£ ^> (2 — £)d. As in many other instances in field theory or statis- 
tical physics, such a non-perturbative formalism is expected to result from a saddle-point 
technique applied to the path integral controlling the statistics of the field. Physically, that 
would correspond to finding some optimal fluctuation responsible for a given structure func- 
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tion. Not any structure function can be found by this approach but only those with N ^> 1, 
related to the PDF tails which are indeed controlled by rare events. This is a general idea 
of the instanton formalism (see, e.g., Coleman, 1977) adapted for turbulence by Falkovich et 
al. (1996). The case in question is so complicated though that an effective analysis (carried 
out by Balkovsky and Lebedev, 1998) is possible only with yet another large parameter, 
(2 — £)d >> 1, which guarantees that the Lagrangian trajectories are almost deterministic. 
The relation between N£ and (2 — £)d is now arbitrary so one is able to describe both the 
perturbative and non-perturbative domains. Unfortunately, a straightforward application 
of this approach to the path integral over the velocity field does not work because of a usual 
problem in saddle-point calculations: the existence of a soft mode makes the integrand non- 
decaying in some direction in the functional space. One ought to integrate over the soft 
mode before the saddle-point approximation is made. Balkovsky and Lebedev identified the 
soft mode as that responsible for the slow variations of the direction of the main stretch- 
ing. Since the structure functions are determined only by the modulus of the distance then 
an effective integration over the soft mode simply corresponds to passing from the velocity 
to the absolute value of the Lagrangian separation as the integration variable in the path 
integral. This can be conveniently done by introducing the scalar variable 

Vl2 = (2 - t)- l d t Rl^ = R^RIM - v\) . (164) 

For the Kraichnan velocity field, rji 2 has the nonzero mean (7712) = —D and the variance 

2D 

«»7l2(tl)W*2))) = -^-912,34 -t 2 ) ■ (165) 

The explicit dependence of the ^12,34 function on the particle distances will not be needed 
here and can be found in the original paper (Balkovsky and Lebedev, 1998). Any average 
over the statistics of the Lagrangian distances can be written in terms of a Martin-Siggia- 
Rose path integral / VRVm exp(tl R ), with the action 



I R = J dt J dridr 2 m 12 



'dtRl^ „\ iD 1 , 
2 _ 1 + D I + -j-dr 3 dr A g 12i34 m 34 
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(166) 



The auxiliary field conjugated to R i2 is denoted by mi 2 = m(t, r±, r 2 ). Note that the second 
(nonlinear) term in (166) vanishes both as (2 — £)d — > oo and £ — > 0. The moments of any 
linear functional of the scalar ■& — J dr(3{r)6(r) are then expressed as 



where T\ = \ J dtdr\ dr 2 xiRu) (3(ri) (3(r 2 ) . To obtain the structure functions, one should 
in principle take for (3 differences of 5 functions. This would however bring diffusive effects 
into the game. To analyze the scaling behavior, it is in fact enough to consider any observable 
where the reducible components in the correlation functions are filtered out. A convenient 
choice is /3(ri) = 5\ (r\ — \ ) — 8\ (ri + |), where the smeared function 5\(r) has a width 
A -1 and satisfies the normalization condition / dr5^(r) = 1. The diffusive effects may be 
disregarded provided the width is taken much larger than r&. 
The saddle-point equations for the integral (167) are 



-iR\ 2 i d t mi2 = — dr 3 dr± * 2 ' 34 [m 12 m u + 2mi 3 m 24 ] + — x\Ri2)P(r 1 )(3(r 2 ) , (169) 



The two boundary conditions are R± 2 (t = 0) = \r 1 — r 2 \ and m 12 — > as t — > —00. The 
variables R\ 2 and m\ 2 are a priori two fields, i.e. they depend on both t and r. In fact, 
the problem can be shown to reduce effectively to two degrees of freedom: describing 
the separation of two points, and R + , describing the spreading of a cloud of size A around a 
single point. It follows from the analysis in (Balkovsky and Lebedev, 1998) that there are two 
different regimes, depending on the order of the moments considered. At iV < (2 — £)d/ (2£) 
the values of R + and i?_ are very close during most of the evolution and different fluid 
particles behave similarly. For higher moments, R + and R- differ substantially throughout 
the evolution. The fact that different groups of fluid particles move in a very different way 




(167) 




(168) 



with the two extremal conditions on the parameters d and y: 




(170) 
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might be interpreted as the signature of the strong fronts in the scalar field that are discussed 
in Section III. F. The final result for the scaling exponents is: 

( N = N(2-0/2-2£N 2 /2d at N < (2 — g)d/(2£) , (171) 
Cat = (2 - 2 d/(80 at N > (2 — £)d/ (2£) . (172) 

These expressions are valid when the fluctuations around instanton give negligible contri- 
bution which requires N ^> 1 and (2 — £)d ^> 1, while the relation between d and N is 
arbitrary. The exponents depend quadratically on the order and then saturate to a con- 
stant. The saturation for the Kraichnan passive scalar model had been previously inferred 
from a qualitative argument by Yakhot (1997) and from an upper bound on ( N by Chertkov 
(1997). The relevance of the phenomenon of saturation for generic scalar turbulence is 
discussed in Section III. F. 

3. Anomalous scaling for magnetic fields 

Magnetic fields transported by a Kraichnan velocity field display anomalous scaling al- 
ready at the level of the second-order correlation functions. For a scalar, (9 2 ) is conserved 
and its flow across the scales fixes the dimensional scaling of the covariance found in Sec- 
tion Ill.C.lii). For a magnetic field, this is not the case. The presence of an anomalous 
scaling for the covariance C^ir^t) = (B' l (r,t) B j (0,t)) becomes quite intuitive from the La- 
grangian standpoint. We have seen in Section II. E that the zero modes are closely related 
to the geometry of the particle configurations. For a scalar field, the single distance involved 
in the two-particle separation explains the absence of anomalies at the second order. The 
magnetic field equation (95) for k — is the same as for a tangent vector, i.e. the separation 
between two infinitesimally close particles. Although C£ involves again two Lagrangian par- 
ticles, each of them is now carrying its own tangent vector. In other words, we are somehow 
dealing with a four-particle problem, where the tangent vectors bring the geometric degrees 
of freedom needed for the appearance of nontrivial zero modes. The compensation mech- 
anism leading to the respective two-particle integral of motion is now due to the interplay 
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between the interparticle distance and the angular correlations of the vectors carried by 
the particles. The attractive feature of the problem is that the anomaly can be calculated 
nonperturbatively. 

Specifically, consider the equation (95) for the solenoidal magnetic field B(r,t) and 
assume the Kraichnan correlation function (48) for the Gaussian incompressible velocity. 
We first analyze the isotropic sector (Vergassola, 1996). The first issue to be addressed 
is the possibility of a stationary state. For this to happen, there should be no dynamo 
effect, i.e. an initial condition should relax to zero in the absence of injection. Let us show 
that this is the case for £ < 1 in 3d. As for a scalar, the ^-correlation of the velocity 
leads to a closed equation for the pair correlation function d t Cl = Mf^Cf. The isotropy 
and the solenoidality of the magnetic field permit to write Cl{r,t) in terms of its trace 
H(r,t) only. This leads to an imaginary time Schrodinger equation for the "wave function" 
•0(r, t) — K+I ^ ir Jq H(p,t) p 2 dp (Kazantsev, 1968). The energy eigenstates ip{r)e~ Et ) into 
which if)(r,t) may be decomposed, satisfy the stationary equation 

fp^l + m(r ) [E _ [/( r )] ^( r ) = o (173) 
dr A 

of a quantum particle of variable mass m{r) in the potential U (r). The presence of a dynamo 
effect is equivalent to the existence of negative energy levels. Since the mass is everywhere 
positive, it is enough to look for bound states in the effective potential V = mil. The 
detailed expressions of the mass and the potential can be found in (Vergassola, 1996). Here, 
it is enough to remark that V(r) is repulsive at small scales and has a quadratic decay in 
the inertial range with a prefactor 2 — 3/2£ — 3/4£ 2 . It is known from quantum mechanics 
textbooks that the threshold for bound states in an attractive potential — c/r 2 is c = 1/4. 
The absence of a dynamo effect for £ < 1 follows immediately from the expression of the 
potential prefactor. 

The stability just described implies that, in the presence of a forcing term in (95), the 
magnetic field covariance will relax to a time-independent expression at long times. We 
may then study the spatial scaling properties at the stationary state. The energy E in 
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(173) should be set to zero and we should add the corresponding forcing term in the first 
equation. Its precise form is not important here as any anomalous scaling is known to 
come from the zero modes of the operator M = d 2 /dr 2 — V(r). The behavior — c/r 2 of the 
potential in the inertial range implies that the operator is scale-invariant and its two zero 
modes behave as power laws with exponents 1/2 1 ± y/\ — Ac . The zero mode with the 
smaller exponent is not acceptable due to a singular behavior in the dissipation range of the 
corresponding correlation function. The remaining zero mode dominates the inertial-range 
behavior H{r) oc r 72 with 



72 = -(3 + 0/2 + 3/2^1 - ^ + 2)/3. (174) 

Note that dimensional arguments based on a constant flux of A 2 (with the vector potential 
defined by B = V x A) would give 72 = — £. This is indeed the result in the limits of small 
£ and large space dimensionality and for the 2d case. For the latter, the vector potential is 
reduced indeed to a scalar whose equation coincides with the advection-diffusion equation 
(92). Note that 72 < — £, that is the zero mode provides for correlation functions that are 
{l/r)~ l2 ~^ times larger than what dimensional arguments would suggest. 

The zero mode dominating the stationary 2-point function of B is preserved by the 
unforced evolution. In other words, if C^ir,^) is taken as the zero mode of M. then the 
correlation function does not change with time. A counterpart to that is the existence of a 
statistical Lagrangian invariant that contains both the distance between the fluid particles 
and the values of the fields: lit) = (B k (R 1 )B l (R 2 )Z^(R 12 )) (Celani and Mazzino, 2000). 
Here is a zero mode of the operator adjoint to M.. The scaling dimension of such zero 
mode is 7 2 + 2 > 0. The appearance of the adjoint operator has a simple physical reason. To 
calculate the correlation functions of the magnetic field, the tangent vectors attached to the 
particles evolve forward in time while, as for the scalar, their trajectories should be traced 
backward in time. Adjoint objects naturally appear when we look for invariants where all 
the quantities run in the same direction of time. The conservation of I(t) is due to the 
power-law growth of Z + with time being offset by the decorrelation between the directions 
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of the B vectors along the separating trajectories. 

Up to now, we have been considering the covariance in the isotropic sector. The scaling 
exponents in the nonisotropic sectors can also be calculated nonperturbatively (Lanotte and 
Mazzino, 1999; Arad, Biferale et al, 2000). The problem is analogous to that for the scalar 
in Section II.E.5, solved by the expression (88). Here, the calculation is more involved 
since in each sector j, m of SO(3) there are nine independent tensors into which C$ may 
be decomposed. Their explicit expression may be found in (Arad et al, 1999). As in the 
isotropic case, it is shown that no dynamo takes place for £ < 1 and the scaling properties 
at the stationary state can then be calculated. For odd j, one has for example: 

ll = -(3 + 0/2 + l/2 v /l-10£ + e + 2j(j + l)(£ + 2). (175) 

The expression for the other sectors can be found in (Lanotte and Mazzino, 1999; Arad, 
Biferale et al, 2000). An important point (we shall come back to it in Section III.F) is that 
the exponents increase with j: the more anisotropic the contribution, the faster it decays 
going toward the small scales. 

Higher-order correlation functions of the magnetic field also obey closed equations. Their 
analysis proceeds along the same lines as for the scalar, with the additional difficulty of the 
tensorial structure. The extra-terms in the equations for the correlation functions of the 
magnetic field come from the presence of the stretching term B ■ Vv. Since the gradient 
of the velocity appears, they are all proportional to £. On the other hand, the stabilizing 
effective-diffusivity terms do not vanish as £ — > 0. We conclude that the forced correlation 
functions at any arbitrary yet finite order will relax to time-independent expressions for 
sufficiently small £. It makes then sense to consider their scaling properties at the stationary 
state, as it was done in (Adzhemyan and Antonov, 1998). For the correlation functions 
((B(r, t) ■ B(0, t)) N ) oc r 7 iv , the perturbative expression in £ reads 7^ = — N£ — 2jV ^ 2 1 ^ + 
0(£ 2 ), demonstrating the intermittency of the magnetic field distribution. 

Let us conclude by discussing the behavior of the magnetic helicity {A ■ B), considered 
in (Borue and Yakhot, 1996). This quantity is conserved in the absence of the molecular 
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diffusivity, as can be easily verified using (95) and the equation for A: 

d t A = vxB -V0+ kV 2 A. (176) 

The function <fi may be fixed by the choice of a specific gauge, e.g. V • A = 0. The 
spatial behavior of the helicity correlation functions (A(r,t) ■ B(0,t)) is derived using (95) 
and (176) and averaging over the Gaussian velocity with the variance (48). The resulting 
equation coincides with that for the scalar covariance (145), implying the dimensional scaling 
r 2 ~£ and a constant helicity flux. We conclude that it is possible to have coexistence of 
normal and anomalous scaling for different components of the correlation tensor of a given 
order. Note also that the helicity correlation functions relax to a stationary form even for 
£ > 1, i.e. when the magnetic correlation functions do not. The increase of the magnetic 
field magnitude is indeed accompanied by a modification of its orientation and the quasi- 
orthogonality between A and B ensures the stationarity of the helicity correlation functions. 
For a helical velocity, considered in (Rogachevskii and Kleeorin, 1999), the magnetic and 
the helicity correlators are coupled via the so-called a-effect, see e.g. (Moffatt, 1978), and 
the system is unstable in the limit n — > considered here. 

D. Lagrangian numerics 

The basic idea of the Lagrangian numerical strategy is to calculate the scalar correlation 
functions using the particle trajectories. The expressions (101) and (134) naturally provide 
for such a Lagrangian Monte-Carlo formulation: the N Lagrangian trajectories are gener- 
ated by integrating the stochastic equations (5), the right hand side of (101) and (134) is 
calculated for a large ensemble of realizations and averaged over it. If we are interested in 
correlation functions of finite order, the Lagrangian procedure involves the integration of 
a few differential equations. This is clearly more convenient than having to deal with the 
partial differential equation for the scalar field. The drawback is that quantities involving a 
large number of particles, such as the tails of the PDF's, are not accessible. Once the corre- 
lation functions have been measured, their appropriate combinations will give the structure 
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functions. For the second-order S 2 (r,t) two different configurations of particles are needed. 
One corresponding to {9 2 }, where the particles are at the same point at time t, and another 
one corresponding to (0(r)0(O)), where they are spaced by r. For the 2n-th order structure 
function, n + 1 particle configurations are needed. 

Another advantage of the Lagrangian method is that it gives direct access to the scaling in 
£ of the structure functions (157), that is to the anomalous dimensions A^ = — (n ■ 

The quantities (101) and (134) for various €s can indeed be calculated along the same 
Lagrangian trajectories. That is more efficient than measuring the scaling in r, i.e. changing 
the final positions of the particles and generating a new ensemble of trajectories. 

1. Numerical method 

The Lagrangian method as presented up to now might be applied to any velocity field. 
The situation with the Kraichnan model is simpler in two respects. First, the velocity 
statistics is time-reversible and the Lagrangian trajectories can be generated forward in 
time. Second, the velocity fields at different times are independent. Only (N — l)d random 
variables are needed at each time step, corresponding to the velocity increments at the 
location of the N particles. The major advantage is that there is no need to generate the 
whole velocity field. Finite-size effects, such as the space periodicity for pseudo-spectral 
methods, are thus avoided. 

The Lagrangian trajectories for the Kraichnan model are conveniently generated as fol- 
lows. The relevant variables are the interparticle separations, e.g. R n N = R n — R^ for 
n = 1, . . . N — 1. Their equations of motion are easily derived from (5) and conveniently 
discretized by the standard Euler-Ito scheme of order 1/2 (Kloeden and Platen, 1992) 

R nN (t + At) - R nN (t) = VAt (V n + V2^W n ) , (177) 

where At is the time step. The quantities V n and W n are ri-dimensional Gaussian, indepen- 
dent random vectors generated at each time step. Both have zero mean and their covariance 
matrices follow directly from the definition (48) of the Kraichnan velocity correlation: 
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(KV£) = d^{R nN ) + d^{R mN ) - d^{R nN - R mN ), (W^Wi) = (1 + S nm )^ . (178) 

The most convenient numerical procedure to generate the two sets of vectors is the classical 
Cholesky decomposition method (Ralston and Rabinowitz, 1978). The covariance matrices 
are triangularized in the form M M T and the lower triangular matrix M is then multiplied 
by a set of (N — l)d Gaussian, independent random variables with zero mean and unit 
variance. The resulting vectors have the appropriate correlations. 

Various possibilities to extract the anomalous scaling exponents are available for the 
Kraichnan model. The straightforward one, used in (Frisch et al, 1998), is to take the 
forcing correlation close to a step function (equal to unity for r/i < 1 and vanishing oth- 
erwise). The correlation functions (134) involve then the products of the average residence 
times of couples of particles at distances smaller than I. An alternative method is based 
on the shape dynamics discussed in Sect. II. E. 4. Measuring first-exit times and not resi- 
dence times gives an obvious advantage in computational time. As stressed in (Gat et al, 
1998), the numerical problem here is to measure reliably the contributions of the irreducible 
zero modes, masked by the fluctuations of the reducible ones. The latter were filtered out 
by Celani et al. (1999) taking various initial conditions and combining appropriately the 
corresponding first-exit times. A relevant combination for the fourth order is for example 
7^(0,0,0,0) - 47^(r , 0,0,0) + 3%(r , r , 0, 0), where % is the first time the size of the 
particle configuration reaches i. 

2. Numerical results 

We shall now present the results for the Kraichnan model obtained by the Lagrangian 
numerical methods just discussed. 

The fourth-order anomaly 2( 2 — d vs the exponent £ of the velocity field is shown in Fig. 5 
for both 2d and 3d (Frisch et al, 1998 and 1999). A few remarks are in order. First, the 
comparison between the 3d curve and the prediction 4£/5 for small £'s provides direct support 
for the perturbation theory discussed in Sect. II. E. 5. Similar numerical support for the 
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expansion in 1/d has been obtained in (Mazzino and Muratore-Ginanneschi, 2001). Second, 
the curve close to £ = 2 is fitted with reasonable accuracy by a(2-£) + &(2-0 3/2 for a = 0.06 
and b = 1.13. That is compatible with an expansion in powers of (2 — £) 1 / 2 (Shraiman and 
Siggia, 1996), where the first term is ruled out by the Schwartz inequality ^4 < 2£ 2 = 2(2— £). 
Third, remark that the anomalies are stronger in 2d than in 3d and their maximum shifts 
towards smaller £ as the dimension decreases. The former remark is in agreement with the 
general idea which emerged in previous Sections that intermittency is associated with the 
particles staying close to each other for times longer than expected. It is indeed physically 
quite sensible that those processes are favored by lowering the space dimensionality. The 
second remark can be qualitatively interpreted as follows (Frisch et al, 1998). Consider 
scalar fluctuations at a given scale. The smaller-scale components of the velocity act like an 
effective diffusivity whilst its larger-scale components affect the scalar as in the Batchelor 
regime. Neither of them leads to any anomalous scaling of the scalar. Those non-local 
interactions are dominant as £ — > and £ — > 2, respectively. For intermediate values of £ 
the velocity components having a scale comparable to that of the scalar fluctuations become 
important and intermittency is produced. The strongest anomalies are attained when the 
relevant interactions are mostly local. To qualitatively explain how the maximum of the 
anomalies moves with the space dimensionality, it is then enough to note that the effective 
diffusivity increases with d but not the large-scale stretching. As for the dependence on 
the order of the moments, the maximum moves toward smaller £ as N increases, see the 
3d curves for the sixth-order anomaly 3(^ — C6 (Mazzino and Muratore-Ginanneschi, 2001). 
It is indeed natural that higher moments are more sensitive to multiplicative effects due to 
large-scale stretching than to additive effects of small-scale eddy diffusivity. 

Let us now discuss the phenomenon of saturation, i.e. the fact that (n tend to a constant 
at large N. The orders where saturation is taking place are expected to increase with £ and 
diverge as £ — > 0. It is then convenient to consider small values of 2 — £. On the other hand, 
approaching the Batchelor limit too closely makes non-local effects important and the range 
of scales needed for reliable measurements becomes huge. A convenient trade-off is that 
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considered in (Celani et al, 2000a) with the 3d cases 2 — £ = 0.125,0.16 and 0.25. For the 
first value of £ it is found there that the fourth and the sixth-order exponents are the same 
within the error bars. The curves for the other £ values show that the order of saturation 
increases with 2 — £, as expected. How do those data constrain the (n curve? It follows from 
the Holder inequalities that the curve for N > 6 must lie below the straight line joining ( 4 
and ^6- Furthermore, from the results in Section III. A. 3 we know that the spatial scaling 
exponents in the forced and the decaying cases are the same and independent of the scalar 
initial conditions. For the unforced equation (92), the maximum value of 9 cannot increase 
with time. Taking an initial condition with a finite maximum value, we have the inequality 
-Sjv < (2max9) p S N ^ p (r,t). We conclude that the ( N curve cannot decrease with the 
order. The presence of error bars makes it, of course, impossible to state rigorously that 
the C/v curves tend to a constant. It is however clear that the combination of the numerical 
data in (Celani et al, 2000a) and the theoretical arguments discussed in Section III.C.2 
leaves little doubt about the saturation effect in the Kraichnan model. The situation with 
an arbitrary velocity field is the subject of Section III.F. 

E. Inverse cascade in the compressible Kraichnan model 

The uniqueness of the Lagrangian trajectories discussed in Section II. D for the strongly 
compressible Kraichnan model has its counterpart in an inverse cascade of the scalar field, 
that is in the appearance of correlations at larger and larger scales. Moreover, the absence 
of dissipative anomaly allows to calculate analytically the statistics of scalar increments and 
to show that intermittency is suppressed in the inverse cascade regime. In other words, the 
scalar increment PDF tends at long times to a scale-invariant form. 

Let us first discuss the simple physical reasons for those results. The absence of a 
dissipative anomaly is an immediate consequence of the expression (101) for the scalar 
correlation functions. If the trajectories are unique, particles that start from the same point 
will remain together throughout the evolution and all the moments (9 N )(t) are preserved. 
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Note that the conservation laws are statistical: the moments are not dynamically conserved 
in every realization, but their average over the velocity ensemble are. 

In the presence of pumping, fluctuations are injected and a flux of scalar variance toward 
the large scales is established. As explained in Section III.B.3, scalar correlation functions 
at very large scales are related to the probability for initially distant particles to come close. 
In strongly compressible flow, the trajectories are typically contracting, the particles tend to 
approach and the distances will reduce to the forcing correlation length I (and smaller) for 
long enough times. Strong correlations at larger and larger scales are therefore established 
as time increases, signaling the inverse cascade process. 

The absence of intermittency is due to the fact that the iV-th order structure function 
is dynamically related to a two-particles process. Correlation functions of the iV-th or- 
der are generally determined by the evolution of iV-particle configurations. However, the 
structure functions involve initial configurations with just two groups of particles separated 
by a distance r. The particles explosively separate in the incompressible case and we are 
immediately back to an iV-particle problem. Conversely, the particles that are initially in 
the same group remain together if the trajectories are unique. The only relevant degrees of 
freedom are then given by the intergroup separation and we are reduced to a two-particle 
dynamics. It is therefore not surprising that the scaling behaviors at the various orders are 
simply related in the inverse cascade regime. 

Specifically, let us consider the equations for the single-point moments (9 N )(t). Since the 
moments are conserved by the advective term, see (101), their behavior in the limit k — > 
(nonsingular now) is the same as for the equation d t 6 = <p. It follows that the single-point 
statistics is Gaussian, with (9 2 ) coinciding with the total injection $(0)t by the forcing. 

The equation for the structure functions Sjy(r,t) is derived from (143). That was im- 
possible in the incompressible case since the diffusive terms could not be expressed in terms 
of Sn (another sign of the dissipative anomaly). No such anomaly exists here so we can 
disregard the diffusion term and simply derive the equation for Sn from (143). This is the 
central technical point allowing for the analytical solution. The equations at the various 
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orders are recast in a compact form via the generating function Z(X ;r,t) = (e tX Ard ) for the 
scalar increments A r 9 = 9(r,t) — 9(0, t). The equation for the generating function is 

d t Z(X ;r,t) = M Z(X ; r, f) + X 2 ($(0) - $(r/£)) Z(X ;r,t), (179) 

where the operator M was defined in (62) and Z = 1 at the initial time. Note that M is 
the restriction of M.2-, signaling the two-particle nature of the dynamics at any order. The 
stationary solution for Z depends on the pumping, but two different regions with a universal 
behavior can be identified. 

i). Large scales r £. In this region, $(r) in (179) can be neglected and the generating 
function is an eigenfunction of M with eigenvalue A 2 $(0). Introducing a new variable 
r (2-£)/2^ the equation (179) is transformed into a Bessel form and solved analytically. The 
solution takes a scale-invariant form Z(\X \ r^ 2 ~^^ 2 ) whose detailed expression may be found 
in (Gaw§dzki and Vergassola, 2000). It follows that S 2 (r) oc r 2 ^. The associated scalar 
flux is calculated as in (148) and turns out to be constant in the upscale direction, the 
footprint of the inverse cascade. The scale-invariant form of the generating function signals 
that no anomalous scaling is present in the inverse cascade regime. As we shall discuss in 
Section IV. B. 2, the phenomenon is not accidental and other physical systems with inverse 
cascades share the same property. Note that, despite its scale-invariance, the statistics of 
the scalar field is strongly non-Gaussian. The expression for the scalar increment PDF 
is obtained by the Fourier transform of the generating function. The tails of the PDF 
decay algebraically with the power —(26 +1), where b = 1 + (7 — d)/(2 — £) and 7 was 
defined in (62). The slow decay of the PDF renders moments of order N > 2b divergent 
(as t N l 2 ~ h r 6 ( 2 ~^) when time increases. A special case is that of smooth velocities £ = 2, 
considered in (Chertkov et ai, 1998). The PDF of scalar increments reduces then to the 
same form as in the direct cascade at small scales. For amplitudes smaller than \n(r/£) the 
PDF is Gaussian. The reason is the same as in Section III.B.2: the time to reach the integral 
scale £ from an initial distance r ^> £ is typically proportional to ln(r/£); the fluctuations of 
the travel time are Gaussian. For larger amplitudes, the PDF has an exponential tail whose 
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exponent depends on the whole hierarchy of the Lyapunov exponents, as in the smooth 
incompressible case. 

ii). Small scales r <C £. Contrary to large scales, the scalar increments are now strongly 
intermittent. The structure functions of integer orders are all dominated by the zero mode of 
the M operator scaling as r b< - 2 ^\ with b defined in the previous paragraph. The exponent 
of the constant flux solution oc r 2 ~^ crosses that of the zero mode at the threshold of 
compressibility b = 1 for the inverse cascade. 

Let us now consider the role of an infrared cutoff in the inverse cascade dynamics. The 
natural motivation is the quasi-stationarity of the statistics: due to the excitation of larger 
and larger scales, some observables do not reach a stationary form. It is therefore of interest 
to analyze the effects of physical processes, such as friction, acting at very large scales. The 
corresponding equation of motion is 

d t 6 + v- V9 + a9 - kV 2 9 = <p, (180) 

and we are interested in the limit a — > 0. For nonsmooth velocities, the friction and the 
advection balance at a scale rj f ~ aT 1 ^ 2- ^, much larger than £ as a — > 0. The smooth case 
£ = 2 is special as no such scale separation is present and it will be considered at the end of 
the Section. For nonsmooth flows, the energy is injected at the integral scale £, transferred 
upwards in the inertial range and finally extracted by friction at the scale r) f . 

The friction term in (180) is taken into account by noting that the field exp(a t) 9 satisfies 
the usual passive scalar equation with a forcing exp(at)f. It follows that the previous 
Lagrangian formulae can be carried over by just introducing the appropriate weights. The 
expression (154) for the iV-point correlation function becomes, for instance, 

C N (r,t) = fds [e-W'ltftr; r>) £ C^r',, . x . ., ; s) *(\r' n - r' m \/£) dr. (181) 

JO J n<m n m 

From (181), one can derive the equations for (9 N )(t) and the structure functions Sjv(r, t) 
and analyze them (Gaw§dzki and Vergassola, 2000). The single-point moments are finite 
and the scalar distribution is Gaussian with (9 2 ) = $(0)/2a. The structure functions of 
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order N < 2b are not affected by friction. The orders that were previously diverging are 
now finite and they all scale as r b ( 2 ~^ in the inertial range, with a logarithmic correction for 
N = 2b. The algebraic tails that existed without friction are replaced by an exponential fall 



off for amplitudes larger than ^/$(0)/a. The saturation of the exponents comes from the 
fact that the moments of order N/2 > b are all dominated by the contribution near the cut. 

Let us conclude by considering the smooth case £ = 2, where the velocity increments 
scale linearly with the distance. The advective term v ■ V has zero dimension, like the 
friction term. As first noted by Chertkov (1998, 1999), this naturally leads to an anomalous 
scaling and intermittency. Let us consider for example the second-order correlation function 
C2(r,t) = (9(r,t)9(0,t)). Its governing equation is derived from (181): 



with the same M operator as in (179). At large scales r ^> £, the forcing term is negligible 
and we look for a stationary solution. Its nontrivial decay C 2 oc (r/l)~ A2 is due to the zero 
mode arising from the balance between the M operator and the friction term: 



where 7 is as in (62) and p is the compressibility degree of the velocity. The notation A 2 is 
meant to stress that dimensional arguments would predict an exponent zero. Higher-order 
connected correlation functions also exhibit anomalous decay laws. Similar mechanisms for 
anomalous scaling and intermittency for the 2c? direct enstrophy cascade in the presence of 
friction are discussed in Section IV.B.l. 

F. Lessons for general scalar turbulence 

The results for spatially nonsmooth flows have mostly been derived within the frame- 
work of the Kraichnan model. Both the forcing and the velocity were Gaussian and short- 
correlated in time. As discussed in Section III.C.l, the conditions on the forcing are not 
crucial and may be easily relaxed. The scaling properties of the scalar correlation functions 




dt C 2 = MC 2 + $(r)C 2 - 2aG 



2- 



(182) 




(183) 
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are universal with respect to the forcing, i.e. independent of its details, while the constant 
prefactors are not. The situation with the velocity field is more interesting and nontrivial. 
Even though a short-correlated flow might in principle be produced by an appropriate forc- 
ing, all the cases of physical interest have a finite correlation time. The very existence of 
closed equations of motion for the particle propagators, which we heavily relied upon, is then 
lost. It is therefore natural to ask about the lessons drawn from the Kraichnan model for 
scalar turbulence in the generic situation of finite-correlated flows. The existing numerical 
evidence is that the basic mechanisms for scalar intermittency are quite robust: anomalous 
scaling is still associated with statistically conserved quantities and the expansion (75) for 
the multiparticle propagator seems to carry over. The specific statistics of the advecting flow 
affects only quantitative details, such as the numerical values of the exponents. The general 
consequences for the universality of the scalar statistics and the decay of the anisotropies 
are presented in what follows. We also show that the phenomenon of saturation, discussed 
in Sections III.C.2 and III.D.2 for the Kraichnan model, is quite general. 

A convenient choice of the velocity v to investigate the previous issues is a 2d flow gener- 
ated by an inverse energy cascade (Kraichnan, 1967). The flow is isotropic, it has a constant 
upscale energy flux and is scale-invariant with exponent 1/3, i.e. without intermittency 
corrections as shown both by experiments and numerical simulations (Paret and Tabeling, 
1998; Smith and Yakhot, 1993; Boffetta et al, 2000). The inverse cascade flow is thus akin 
to the Kraichnan ensemble, but its correlation times are finite. 

Let us first discuss the preserved Lagrangian structures. The simplest nontriv- 
ial case to analyze anomalous scaling is the third-order correlation function C 3 (r) = 
(0(ri, i) 9(r 2 ,t) #(r 3 ,t)). The function is non-zero only if the symmetry 9 —9 is broken, 
which often happens in real systems via the presence of a mean scalar gradient (9) = g ■ r. 
The function C 3 (r) depends then on the size, the orientation with respect to g and the 
shape of the triangle defined by t*i, r 2 and r 3 . For a scalar field advected by the 2d in- 
verse energy cascade flow, the dependence on the size R of the triangle is a power law with 
anomalous exponent £3 = 1.25, smaller than the dimensional prediction 5/3 (Celani et al, 
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2000a). To look for statistical invariants under the Lagrangian dynamics, let us take a 
translation-invariant function f(r) of the N points r n and define its Lagrangian average as 
in (72), i.e. as the average of the function calculated along the Lagrangian trajectories. In 
the 2d inverse energy cascade, the distances grow as \t\ 3 ^ 2 and the Lagrangian average of a 
scaling function of positive degree a is expected to grow as |t| 3cr//2 . The numerical evidence 
presented in (Celani and Vergassola, 2001) is that the anomalous scaling is again due to 
statistical integrals of motion: the Lagrangian average of the anomalous part of the corre- 
lation functions remains constant in time. The shape of the figures identified by the tracer 
particles plays again a crucial role: the growth oc |t| 3< > 3 / 2 of the size factor R 1 ^ 3 in C$(r) is 
compensated by its shape dependence. As we indicated in Section II.E.4, anomalous scaling 
reflects slowed-down separations among subgroups of particles and the fact that triangles 
with large aspect ratios live much longer than expected. It is also immediate to provide 
an example of slow mode, see (75), for the two particle dynamics. The Lagrangian average 
of (g ■ 7*12) is obviously preserved as its time derivative is proportional to ((vi — v 2 )) = 0. 
Its first slow mode is given by (g ■ 7*12) r^ 3 , whose Lagrangian average is found to grow as 
\t\ (that is much slower than |t| 5 ^ 2 ) at large times. In the presence of a finite volume and 
boundaries, the statistical conservation laws hold as intermediate asymptotic behaviors, as 
explicitly shown in (Arad et al, 2001). 

An important consequence is about the decay rate of the anisotropies. As already men- 
tioned, isotropy is usually broken by the large-scale injection mechanisms. One of the 
assumptions in the Oboukhov-Corrsin reformulation of the Kolmogorov 1941 theory for the 
passive scalar is that the statistical isotropy of the scalar is restored at sufficiently small 
scales. Experiments do not confirm those expectations. Consider, for example, the case 
where a mean scalar gradient g is present. A quantitative measure of the degree of anisotropy 
is provided by odd-order structure functions or by odd-order moments of g ■ V6. All these 
quantities are identically zero for isotropic fields. The predictions of the Oboukhov-Corrsin 
theory for the anisotropic situations are the following. The hyperskewness of the scalar 
increments S2 n +i(r) / S% +1 ^ 2 (r) should decay with the separation as r 2 / 3 . The corresponding 
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behavior of the scalar gradient hyperskewness with respect to the Peclet number should be 
Pe -1 / 2 . In fact, the previous quantities are experimentally found to remain constant or even 
to increase with the relevant parameter (Gibson et al, 1977; Mestayer, 1982; Sreenivasan, 
1991; Mydlarski and Warhaft, 1998). There is therefore no restoration of isotropy in the 
original Kolmogorov sense and the issue of the role of anisotropies in the small-scale scalar 
statistics is naturally raised (Sreenivasan, 1991; Warhaft, 2000). The analysis of the same 
problem in the Kraichnan model is illuminating and permits to clarify the issue in terms of 
zero medes and their scaling exponents. For isotropic velocity fields, the correlation func- 
tions may be decomposed according to their angular momentum j, as in Section II. E. 5 Each 
of those contributions is characterized by a scaling exponent ( 3 N . The general expectation, 
confirmed in all the situations where the explicit calculations could be performed, is that 
Cif° > C^ = ° an d that the exponents increase with j. As their degree of anisotropy increases, 
the contributions are less and less relevant at small scales. Note that, in the presence of 
intermittency, the inequality < yC 2 i s s ^ m possible and anisotropies might then have 
dramatic effects for observables whose isotropic dominant contribution is vanishing, such as 
5 , 2ri+i/5 , 2 +1//2 - Rather than tending to zero, they may well increase while going toward the 
small scales, blatantly violating the restoration of isotropy in the original Kolmogorov sense. 
Remark that no violation of the hierarchy in j is implied though. In other words, the degree 
of anisotropy of every given moment does not increase as the scale decreases; if however one 
measures odd moments in terms of the appropriate power of the second one (as is customary 
in phenomenological approaches) then the degree of anisotropy may grow downscales. The 
previous arguments are quite general and compatible with all the existing experimental and 
numerical data for passive scalar turbulence. For Navier-Stokes turbulence, the use of the 
rotation symmetries and the existence of a hierarchy among the anisotropic exponents were 
put forward and exploited in (Arad et al. 1998, 1999). Numerical evidence for persistence 
of anisotropies analogous to those of the scalar fields was first presented in (Pumir and 
Shraiman, 1995). 

Let us now discuss the phenomenon of saturation. A snapshot of the scalar field advected 
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by the 2d inverse cascade flow is shown in Fig. 6. A clear feature is that strong scalar 
gradients tend to concentrate in sharp fronts separated by large regions where the variations 
are weak. The scalar jumps across the fronts are of the order of 9 rms = \J (9 2 ), i.e. comparable 
to the largest values of the field itself. Furthermore, the minimal width of the fronts reduces 
with the dissipation scale, pointing to their quasi-discontinuous nature. If the probability 
of having such 6 rms jumps across a separation r decreases as r^°°, then phenomenological 
arguments of the multifractal type suggest a saturation to the asymptotic value £00, see 
(Frisch, 1995). The presence of fronts in scalar turbulence is a very well established fact, 
both in experiments (Gibson et al, 1977; Mestayer, 1982; Sreenivasan, 1991; Mydlarski and 
Warhaft, 1998) and in numerical simulations (Holzer and Siggia, 1994; Pumir, 1994). It is 
shown in the latter work that fronts are formed in the hyperbolic regions of the flow, where 
distant particles are brought close to each other. The other important remark is that fronts 
appear also in the Kraichnan model (Fairhall et al, 1997; Vergassola and Mazzino, 1997; 
Chen and Kraichnan, 1998), despite the ^-correlation of the velocity. What matters for 
bringing distant particles close to each other are indeed the effects cumulated in time. The 
integral of a 5-correlated random process behaves as a Brownian motion in time, whose sign 
is known to have strong persistence properties (Feller, 1950). Even a 5-correlated flow might 
then efficiently compress the particles (locally) and this naturally explains how fronts may 
be formed in the Kraichnan model. It is also clear from the previous arguments that a finite 
correlation time favors the formation of fronts and that the Kraichnan model is somehow the 
most unfavorable case in this respect. The fact that fronts still form and that the saturation 
takes place points to generality of the phenomenon for scalar turbulence. The order of the 
moments and the value where the (n curve flattens out might depend on the statistics of 
the advecting velocity, but the saturation itself should generally hold. Direct evidence for 
the advection by a 2d inverse cascade flow is provided in (Celani et al, 2000a,b). Saturation 
is equivalent to the scalar increment PDF taking the form V(A r 6) = r^°°q(A r 9/9 rms ) for 
amplitudes larger than 6 rms . The tails at various r can thus be all collapsed by plotting 
r -Coo-p ) as shown in Fig. 7. Note finally that the saturation exponent Coo coincides with the 
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fractal codimension of the fronts, see (Celani et al, 2001) for a more detailed discussion. 

As far as the compressible Kraichnan model is concerned, applying even qualitative pre- 
dictions requires much more care than in the incompressible case. Indeed, the compressibility 
of the flow makes the sum of the Lyapunov exponents nonzero and leads to the permanent 
growth of density perturbations described in Sect. III. A. 4. In a real fluid, such growth is 
stopped by the back-reaction of the density on the velocity, providing for a long-time mem- 
ory of the divergence V • v of the velocity along the Lagrangian trajectory. This shows that 
some characteristics of the Lagrangian velocity may be considered short-correlated (like the 
off-diagonal components of the strain tensor), while other are long correlated (like the trace 
of the strain). 

In summary, the situation with the Kraichnan model and general passive scalar turbu- 
lence is much like the motto at the beginning of the review. The interest was originally 
stirred by the closed equations of motion for the correlation functions and the possibility 
of deriving an explicit formula for the anomalous scaling exponents, that are quite specific 
features. Actually, the model turned out to be much richer and capable of capturing the 
basic properties of the Lagrangian tracer dynamics in generic turbulent flow. The major 
lessons drawn from the model, such as the statistical integrals of motion, the geometry 
of the particle configurations, the dynamics in smooth flow, the importance of multipoint 
correlations, the persistence of anisotropies, all seem to have quite general validity. 

IV. BURGERS AND NAVIER-STOKES EQUATIONS 

All the previous sections were written under the assumption that the velocity statistics 
(whatever it is) is given. In this Chapter, we shall describe what one can learn about the 
statistics of the velocity field itself by considering it in the Lagrangian frame. We start from 
the simplest case of Burgers turbulence whose inviscid version describes a free propagation 
of fluid particles, while viscosity provides for a local interaction. We then consider an in- 
compressible turbulence where the pressure field provides for a nonlocal interaction between 
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infinitely many particles. 



A. Burgers turbulence 

The d- dimensional Burgers equation (Burgers, 1974) is a pressureless version of the 
Navier-Stokes equation (1): 

dtV + v-Vv- vV 2 v = f (184) 

for irrotational (potential) velocity v(r, t) and force f(r, t). It is used to describe a variety of 
physical situations from the evolution of dislocations in solids to the formation of large scale 
structures in the universe, see e.g. (Krug and Spohn, 1992; Shandarin and Zel'dovich, 1989). 
Involving a compressible v, it allows for a meaningful (and nontrivial) one-dimensional case 
that describes small-amplitude perturbations of velocity, density or pressure depending on 
a single spatial coordinate in the frame moving with sound velocity, see e.g. (Landau and 
Lifshitz, 1959). Without force, the evolution described by (184) conserves total momentum 
/ v dr. By the substitutions v = V0 and / = Vg the Burgers equation is related to the 
KPZ equation 

<9 t 0+i(V0) 2 -z/V 2 = s. (185) 

governing an interface growth (Kardar et al, 1986). 

Already the one-dimensional case of (184) illustrates the themes that we discussed in 
previous Chapters: turbulent cascade, Lagrangian statistics and anomalous scaling. Under 
the action of a large-scale forcing (or in free decay of large-scale initial data) a cascade 
of kinetic energy towards the small scales takes place. The nonlinear term provides for 
steepening of negative gradients and the viscous term causes energy dissipation in the fronts 
that appear this way. In the limit of vanishing viscosity, the energy dissipation stays finite 
due to the appearance of velocity discontinuities called shocks. The Lagrangian statistics 
is peculiar in such an extremely nonsmooth flow and can be closely analyzed even though 
it does not correspond to a Markov process. Forward and backward Lagrangian statistics 
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are different, as it has to be in an irreversible flow. Lagrangian trajectories stick to the 
shocks. That provides for a strong interaction between the particles and results in an 
extreme anomalous scaling of the velocity field. A tracer field passively advected by such a 
flow undergoes an inverse cascade. 

At vanishing viscosity, the Burgers equation may be considered as describing a gas of 
particles moving in a force field. Indeed, in the Lagrangian frame defined for a regular 
velocity by R = v(R,t), relation (184) becomes the equation of motion of non-interacting 
unit-mass particles whose acceleration is determined by the force: 

R=f(R,t). (186) 

In order to find the Lagrangian trajectory R(t; r) passing at time zero through r it is then 
enough to solve the second order equation (186) with the initial conditions R(0) = r and 
R(0) = v(r, 0). For sufficiently short times such trajectories do not cross and the Lagrangian 
map r i— > R(t; r) is invertible. One may then reconstruct v at time t from the relation 
v(R(t),t) = R(t). The velocity stays potential if the force is potential. At longer times, 
however, the particles collide creating velocity discontinuities, i.e. shocks. The nature and 
the dynamics of the shocks may be understood by treating the inviscid equation as the limit 
of the viscous one. Positive viscosity removes the singularities. As is well known, the KPZ 
equation (185) may be linearized by the Hopf-Cole substitution Z = exp[— (f>/2v\ that gives 
rise to the heat equation in an external potential: 

[d t -vV 2 + lg]Z = 0. (187) 

The solution of the initial value problem for the latter may be written as the Feynman-Kac 
path-integral 



?(r,f)= / exp 



Z(R(0),0) VR. (188) 



R(t)=r 



* r • 2 i 

with the classical action S(R) = J \R +g{r,R) dr of a path R(t) . The limit of vanishing 

o L J 

viscosity selects the least-action path: 
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<j>(r,t) = min \S(R) + <f>(R(0),0) \ . (189) 

R : R(t)=T ' 

Equating to zero the variation of the minimized expression, one infers that the minimizing 
path R(t) = R(r;r,t) is a solution of (186) such that R(0) = v(R(0),0). Taking the 
gradient of (189), one also infers that v(r, t) = R(t). The above procedure extends the short- 
time construction of the solutions of the inviscid Burgurs equation to all times at the cost 
of admitting shocks where v is discontinuous. The velocity still evolves along the solutions 
of (186) with R(0) given by the initial velocity field, but if there are many such solutions 
reaching the same point at the same moment, the ones that do not realize the minimum 
(189) should be disregarded. The shocks arise when there are several minimizing paths. 
At fixed time, shocks are, generically, hypersurfaces which may intersect, have boundaries, 
corners, etc. (Vergassola et al, 1994; Frisch et al, 1999; Frisch and Bee, 2000). This may 
be best visualized in the case without forcing where the equation (189) takes the form 

0(r,t) = min jl(r-r') 2 + <f>(r', 0) } (190) 

with a clear geometric interpretation: <f>(r, t) is the height C of the inverted paraboloid 
C — ^(r — r') 2 touching but not overpassing the graph of 0(r',O). The points of contact 
between the paraboloid and the graph correspond to r' in (190) on which the minimum is 
attained. Shock loci are composed of those r for which there are several such contact points. 

The simplest situation occurs in one dimension with r = x. Here, shocks are located at 
points Xi where the velocity jumps, i.e. where the two-sided limits v ± (x, t) = lim v(x ± e, t) 

e— >0 

are different. The limits correspond to the velocities of two minimizing paths that do not 
cross except at Xj. The shock height Sj = v + {x i) t) — v~(xi,t) has to be negative. Once 
created shocks never disappear but they may merge so that they form a tree branching 
backward in time. For no forcing, the shock positions x« correspond to the inverse parabolas 
C — ^(xi — x') 2 that touch the graph of 0(x',O) in (at least) two points x'^, x 1 ^ > x'~, 
such that v ± (x i ,t) = Xl ~* 1 , see Fig. 8. Below, we shall limit our discussion to the one- 
dimensional case that was most extensively studied in the literature, see (Burgers, 1974; 
Woyczyhski, 1998; E and Vanden Eijnden, 2000a). In particular, the asymptotic long-time 
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large-distance behavior of freely decaying initial data with Gaussian finitely correlated ve- 
locities or velocity potentials has been intensively studied. In the first case, the asymptotic 
solution for x 3 / 2 t~ 1 = 0(1) has the form (190) with (f>(x', 0) representing a Brownian mo- 
tion. In the second case, the solution settles under the diffusive scaling with a logarithmic 
correction x 2 ^ 1 In 1 / 4 x = 0(1) to the form 

<f>(x,t) = nun {^(x-%) 2 + 0,} , (191) 

where (yj,(f>j) is the Poisson point process with intensity e^ dydcj). In both cases explicit 
calculations of the velocity statistics have been possible, see, respectively, (Burgers, 1974; 
Frachebourg and Martin, 2000) and (Kida, 1979; Woyczyhski, 1998). Other asymptotic 
regimes of decaying Gaussian initial data were analyzed in (Gurbatov et ai, 1997). 

The equation of motion of the one-dimensional shocks Xi(t) is easy to obtain even in 
the presence of forcing. To this aim, note that along the shock there are two minimizing 
solutions determining the same function <p and that 

^4>{x h t) = d t (j>{x u t) +x i d x <j)(x i ,t) = ~v ± (x i ,t) 2 + g(xi,t) + i i w ± (x i ,t). (192) 

Equating both expressions, we infer that Xi = ^[v + (xi, t) +v~(xi, £)] = v(xi, t) = Vi , i.e. that 
the shock speed is the mean of the velocities on both sides of the shock. The crucial 
question for the Lagrangian description of the Burgers velocities is what happens with the 
fluid particles after they reach shocks where their equation of motion x = v(x,t) becomes 
ambiguous. The question may be easily answered by considering the inviscid case as a limit 
of the viscous one where shocks become steep fronts with large negative velocity gradients. It 
is easy to see that the Lagrangian particles are trapped within such fronts and keep moving 
with them. We should then define the inviscid Lagrangian trajectories as solutions of the 
equation x = v(x,t), with x understood as the right derivative. Indeed, the solutions of 
that equation clearly move with the shocks after reaching them. In other words, the two 
particles arriving at the shock from the right and the left at a given moment aggregate upon 
the collision. Momentum is conserved so that their velocity after the collision is the mean 
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of the incoming ones (recall that the particles have unit mass) and is equal to the velocity 
of the particles moving with the shock that have been absorbed at earlier times. Note that 
in the presence of shocks the Lagrangian map becomes many-to-one, compressing whole 
space-intervals into the shock locations. 

It is not difficult to write field evolution equations that take into account the presence of 
shocks (Vol'pert, 1967; Bernard and Gaw§dzki, 1998; E and Vanden Eijnden, 2000a). We 
shall do it for local functions of the velocity of the form e Xv ^ x '^ . For positive viscosity, these 
functions obey the equation of motion 

(9 t + A,9 A A- 1 ^-A/)e At ' = -A 2 e(A) , (193) 

where e(A) = u[(d x v ) 2 — A~ 2 <9 2 ] e Xv is the contribution of the viscous term in (184). In the 
limit v — > the dissipation becomes concentrated within the shocks. Using the representa- 
tion e Xv< ^ x ' t ) = e *v + (x,t) q^ x _ Xi (t}} _|_ e Xv fat) g(xi(t) — x) around the shocks, it is easy to 
check that (193) still holds for v = with 

e(X;x,t) = J2F(vi,Si)S(x-Xi(t)), (194) 

i 

concentrated at shock locations. The "form-factor" F(v, s) = —2e Xv A -1 d\ A -1 sinh 4^. 

When the forcing and/or initial data are random, the equation of motion (193) induces 
the Hopf evolution equations for the correlation functions. For example, in the stationary 
homogeneous state (/ e Xv ) = A (e(A)). Upon expanding into powers of A, this relates the 
single-point expectations (v n ) to the shock statistics. The first of these relations says that 
the average force should vanish and the second gives the energy balance (fv) = e v , where 
the mean energy dissipation rate e v = (e(0)) = — p(s 3 )/12 and p = (J2$( x — x i(t))) is 

i 

the mean shock density. Similarly, for the generating function of the velocity increments 
(exp(A Av)) with Av — v± — v 2 , one obtains 

\ 2 d x \- 2 d l (e XAv ) - \(Afe XAv ) = -X 2 (e(X) ie - Xv2 +e Xvi e(-X) 2 ). (195) 

For a Gaussian force with (f(x, t) /(0, 0)) = 5(t)x{%/L), where L is the injection scale, 
(Afe XAv ) = [x(0) — x(Ax/L)] (e XAv ). For a spatially smooth force, the second term 
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in (195) tends to zero when the separation Ax (taken positive) shrinks. In the same 
limit, the quantity inside the expectation value on the right hand side tends to a local 
operator concentrated on shocks, as in (194), but with the height- dependent form-factor 
_ 4e As / 2 A -1 dAA -1 sinh 4^ = — d\\~ 2 (e Xs — 1 — As). Comparing the terms, one infers that 
for = 3,4..., 

lim di ((Av) N ) = p(s N ). (196) 

The first of these equalities is the one-dimensional version of the Kolmogorov flux relation 
(2). The higher ones express the fact that, for small Ax, the higher moments of velocity 
increments are dominated by the contribution from a single shock of height s occurring with 
probability p Ax. That implies the anomalous scaling of the velocity structure functions 

S N (Ax) = ((Av) N ) = p(s N )Ax + o(Ax), (197) 

that is a signature of an extreme intermittency. In fact, due to the shock contributions, all 
moments of velocity increments of order p > 1 scale with exponents ( p = 1. In contrast, for 
< p < 1, ( p = p since the fractional moments are dominated by the regular contributions 
to velocity. Indeed, denoting by £ the regular part of the velocity gradient, one obtains (E 
et al, 1997; E, Khanin et al. 2000; E and Vanden Eijnden, 2000a) 

((\Av\Y) = m p ) (M P + o((AxY) ■ (198) 

The shock contribution proportional to Ax is subleading in that case. 

The Lagrangian interpretation of these results can be based on the fact that the velocity 

is a Lagrangian invariant of the unforced inviscid system. In the presence of the force, 

t 

v(x,t) = v(x(0),0) + J f(x(s),s)ds, (199) 

o 

along the Lagrangian trajectories. The velocity is an active scalar and the Lagrangian 
trajectories are evidently dependent on the force that drives the velocity. One cannot write 
a formula like (135) obtained by two independent averages over the force and over the 
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trajectories. Nevertheless, the main contribution to the distance-dependent part of the 
2-point function (v(x, t) v(x', t)) is due, for small distances, to realizations with a shock 
in between the particles. It is insensitive to a large-scale force and hence approximately 
proportional to the time that the two Lagrangian trajectories ending at x and x' take to 
separate backwards to the injection scale L. With a shock in between x and x' at time t, 
the initial backward separation is linear so that the second order structure function becomes 
proportional to Ax, in accordance with (197). Other structure functions may be analyzed 
similarly and give the same linear dependence on the distance (all terms involve at most two 
trajectories). This mechanism of the anomalous scaling is similar to that of the collinear 
case in Sect. III.B.3. 

Following numerical observations of Chekhlov and Yakhot (1995, 1996), a considerable 
effort has been invested to understand the shape of the PDF's V(Av) of the velocity incre- 
ments and V(£) of the regular part of d x v . The stationary expectation of the exponential 
of £ satisfies for a white-in-time forcing the relation 

[\dl-d x -DX 2 ) (e x Z) = ( Px ) (200) 

that may established the same way as (195). Here D — — \ x"(0) L~ 2 and p x is an operator 
supported on shocks with the form-factor |(e A ^ + +e A ^ ). For the PDF of £ this gives the 
identity (E and Vanden Eijnden, 1999, 2000a): 

(DSf + ^ + 3$) V(0 + / ^e-^( Px ) = 0. (201) 

Various closures for this equation or for (195) have been proposed (Polyakov, 1995; Bouchaud 
and Mezard, 1996; Boldyrev, 1997; Gotoh and Kraichnan, 1998). They all give the right tail 
oc e~£ ^ 3D ^ of the distribution, as determined by the first two terms with derivatives in (201), 
with the power-law prefactors depending on the details of the closure. That right tail was first 
obtained by Feigel'man (1980) and also reproduced for V(Av) by an instanton calculation 
in (Gurarie and Migdal, 1996). The instanton appears to be a solution of a deterministic 
Burgers equation with a force linear in space, see also (Balkovsky et ai, 1997b). The right 
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tail may be also understood from the stochastic equation along the trajectory -§^Ax = Af. 
For Ax much smaller than the injection scale L, the force increment may be linearized 
Af = a(t) Ax. In the first approximation, a(t) is a white noise and we obtain a problem 
familiar from the one-dimensional localization in a 5-correlated potential. In particular, 
z = j- t (Ax) J Ax satisfies z = a(t) — z 2 . This is a Langevin dynamics in the unbounded 
from below potential oc z 3 (Bouchaud and Mezard, 1996). Upon conditioning against escape 
to z — — oo, one gets a stationary distribution for z with the right tail oc e~ z3 /( 3D ) and the 
left tail oc z' 2 . In reality, due to the shock creation, Af / Ax m d x f(x(t),t) is a white 
noise only if we fix the value of v at the end point of the Lagrangian trajectory. This 
introduces subtle correlations which effect the power-law factors in V(£), in particular its 
left tail oc |£|~ Q . Large negative values of £ appear in the vicinity of preshocks, as stressed 
in (E et al, 1997) where the value a = 7/2 was argued for, based on a geometric analysis of 
the preshocks, the birth points of the shock discontinuities. E and Vanden Eijnden (2000a), 
have proved by analysis of realizability conditions for the solutions of (201) that a > 3 
and made a strong case for a = 7/2, assuming that shocks are born with vanishing heights 
and that preshocks do not accumulate. The exponent a = 7/2 was subsequently found in 
the decaying case with random large scale initial conditions both in Id (Bee and Frisch, 
2000) and in higher dimensions (Frisch et al, 2001), and in the forced case when the forcing 
consists of deterministic large scale kicks repeated periodically (Bee et al, 2000). In those 
cases, the statistics of the shocks and there creation process are easier to control than for 
the 5-correlated forcing. The numerical analysis of the latter case clearly confirms, however, 
the prediction a = 7/2. As to the PDF V(Av) of the velocity increment in the decaying 
case and, possibly, in the forced case, it exhibits a crossover from the behavior characteristic 
for the velocity gradients to the one reproducing the behaviors (197) and (198). Note that 
the single-point velocity PDF also has cubic tails \nV(v) oc — \v | 3 (Avellaneda et al, 1995; 
Balkovsky et al, 1997b). The same is true for white-driven Navier-Stokes equation and may 
be generalized for a non-Gaussian force (Falkovich and Lebedev, 1997). 
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The Lagrangian picture of the Burgers velocities allows for a simple analysis of advection 
of scalar quantities carried by the flow. In the inviscid and diffusionless limit, the advected 
tracer satisfies the evolution equation 

d t 9 + v d x 9 = ip, (202) 

where ip represents an external source. As usual, the solution of the initial value problem 
is given in terms of the PDF p(x, t; y, | v) to find the backward Lagrangian trajectory at 
y at time 0, given that at later time t it passed by x, see (100). Except for the discrete 
set of time t shock locations, the backward trajectories are uniquely determined by x. As a 
result, a smooth initial scalar will develop discontinuities at shock locations but no stronger 
singularities. Since a given set of points (x 1 , . . . ,x N ) = x avoids the shocks with probability 
1, the joint backward PDF's of iV trajectories V N (x;y; —t), see (65), should be regular for 
distinct x n and should possess the collapse property (67). This leads to the conservation of 
(9 2 ) in the absence of scalar sources and to the linear pumping of the scalar variance into 
a soft mode when a stationary source is present. Such behavior corresponds to an inverse 
cascade of the passive scalar, as in the strongly compressible phase of the Kraichnan model 
discussed in Section (III.E). 

The Burgers velocity itself and all its powers constitute an example of advected scalars. 
Indeed, the equation of motion (193) may be also rewritten as 

(d t + vd x -Xf)e Xv = (203) 

from which the relation (202) for = v n and ip = nfv™"^ 1 follows. Of course, v n are active 
scalars so that in the random case their initial data, the source terms, and the Lagrangian 
trajectories are not independent, contrary to the case of passive scalars. That correlation 
makes the unlimited growth of (v 2 ) impossible: the larger the value of local velocity, the 
faster it creates a shock and dissipates the energy. The difference between active and passive 
tracers is thus substantial enough to switch the direction of the energy cascade from inverse 
for the passive scalar to direct for the velocity. 
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As usual in compressible flows, the advected density n satisfies the continuity equation 

d t n + d x (vn) = ip (204) 

different from (202) for the tracer. The solution of the initial value problem is given by the 
forward Lagrangian PDF p(y, 0;x,t\v), see (100). Since the trajectories collapse, a smooth 
initial density will become singular under the evolution, with 5-function contributions con- 
centrating all the mass from the regions compressed to shocks by the Lagrangian flow. Since 
the trajectories are determined by the initial point y, the joint forward PDF's V N (y;x;t) 
should have the collapse property (67) but they will also have contact terms in x n 's when 
the initial points y n are distinct. Such terms signal a finite probability of the trajectories 
to aggregate in the forward evolution, the phenomenon that we have already met in the 
strongly compressible Kraichnan model discussed in Sect. (II. D). The velocity gradient d x v 
is an example of an (active) density satisfying equation (204) with ip = d x f. 

In (Bernard and Bauer, 1999), the behavior of the Lagrangian PDF's and the advected 
scalars summarized above have been established by a direct calculation in freely decaying 
Burgers velocities with random Gaussian finitely-correlated initial potentials <fi. 

B. Incompressible turbulence from a Lagrangian viewpoint 

As we learned from the study of passive fields, treating the dissipation is rather easy 
as it is a linear mechanism. The main difficulty resides in proper understanding advection. 
For incompressible turbulence, the problem is even more complicated than for the Burgers 
equation due to spatial nonlocality of the pressure term. The Euler equation may indeed be 
written as the equation 

R = f(R,t)-VP, (205) 

for the Lagrangian trajectories R(t; r) where / is the external force and the pressure field 
is determined by the incompressibility condition V 2 P = —V ■ [v ■ Vw] with v = R and 
the spatial derivatives taken with respect to R. The inversion of the Laplace operator in 
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the previous relation brings in nonlocality via the kernel decaying as a power law. We thus 
have a system of infinitely many particles interacting strongly and nonlocally. In such a 
situation, any attempt at an analytic description looks unavoidably dependent on possible 
simplifications in limiting cases. The natural parameter to exploit for the incompressible 
Euler equation is the space dimensionality, varying between two and infinity. The two- 
dimensional case indeed presents important simplifications since the vorticity is a scalar 
Lagrangian invariant of the inviscid dynamics, as we shall discuss hereafter. The opposite 
limit of the infinite-dimensional Euler equation is very tempting for some kind of mean-field 
approximation to the interaction among the fluid particles but nobody has derived it yet. 
The level of discussion in this section is thus naturally different from the rest of the review: 
as a consistent theory is absent, we present a set of particular arguments and remarks that, 
on one hand, make contact with the previously discussed subjects, and, on other hand, may 
inspire further progress. 

1. Enstrophy cascade in two dimensions 

The Euler equation in any even- dimensional space has an infinite set of integrals of motion 
besides energy. One may indeed show that the determinant of the matrix uj^ = V jV l — Vjf- 7 
is the nonnegative density of an integral of motion, i.e. / F(detcu)dr is conserved for any 
function F. The quadratic invariant /(det u) 2 dr is called enstrophy. In the presence of an 
external pumping injecting energy and enstrophy, it is clear that both quantities may flow 
throughout the scales. If both cascades are present, they cannot go in the same direction: 
the different dependence of energy and enstrophy on the scale prevents their fluxes to be 
both constant in the same interval. A finite energy dissipation would imply an infinite 
enstrophy dissipation in the limit v — > 0. The natural conclusion is that, given a single 
pumping at some intermediate scale and assuming the presence of two cascades, the energy 
and the enstrophy flow toward the large and the small scales, respectively (Kraichnan, 1967; 
Batchelor, 1969). This is indeed the case for the two-dimensional case. 
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In this Section, we shall focus on the 2d direct enstrophy cascade. The basic knowledge 
of the Lagrangian dynamics presented in Sections II.B.l and III.B is essentially everything 
needed. Vorticity in 2d is a scalar and the analogy between vorticity and passive scalar was 
noticed by Batchelor and Kraichnan already in the sixties. Vorticity is not passive though 
and such analogies may be very misleading, as it is the case for vorticity and magnetic field 
in 3d and for velocity and passive scalar in Id Burgers. The basic flux relation for the 
enstrophy cascade is analogous to (132): 

((vi ■ Vi + v 2 ■ V2V1CJ2) = (y?i^2 + ^2^1) = P2 • (206) 

The subscripts indicate the spatial points r x and r 2 and the pumping is assumed to be 
Gaussian with (<p(r, £)0(O, 0)) = 5(t)$(r / L) decaying rapidly for r > L. The constant 
P 2 = $(0), having dimensionality time -3 , is the input rate of the enstrophy uj 1 . Equation 
(206) states that the enstrophy flux is constant in the inertial range, that is for r\ 2 much 
smaller than L and much larger than the viscous scale. A simple power counting suggests 
that the velocity differences and the vorticity scale as the first and the zeroth power of r 12 , 
respectively. That fits the idea of a scalar cascade in a spatially smooth velocity: scalar 
correlation functions are indeed logarithmic in that case, as it was discussed in Sect. III.B. 

Even though one can imagine hypothetical power-law vorticity spectra (Saffman, 1971; 
Moffatt, 1986; Polyakov, 1993), one can argue that they are structurally unstable (Falkovich 
and Lebedev, 1994). Indeed, imagine for a moment that the pumping at L produces the 
spectrum {{uj\ — uj 2 ) p ) oc r^ 2 at r 12 <C L. Regularity of the Euler equation in 2d requires 
( p > 0, see, e.g., (Eyink, 1995) and references therein. In the spirit of the stability theory 
of Kolmogorov spectra (Zakharov et ai, 1992), let us add an infinitesimal pumping at 
some £ in the inertial interval producing a small yet nonzero flux of enstrophy. Small 
perturbations 5u obey the equation d t Su + (v'V)Slj + (5vV)uj = uV 2 5uj. Here, 5v is the 
velocity perturbation related to 5oj. The perturbation 5u has the typical scale £ while 
vorticity u, associated with the main spectrum, is concentrated at L when ( 2 > 0. The 
third term can be neglected as it is (£/L) 2 times smaller than the second one. Therefore, 
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5uj behaves as a passive scalar convected by main turbulence, i.e. the Batchelor regime 
from Section III.B takes place. The correlation functions of the scalar are logarithmic in 
this case for any velocity statistics. The perturbation in any vorticity structure function 
thus decreases downscales slower than the contribution of the main flow. That means that 
any hypothetical power-law spectrum is structurally unstable with respect to the pumping 
variation. Stability analysis cannot of course describe the spectrum downscales where the 
perturbations is getting comparable to the main flow. It is logical though to assume that 
since only the logarithmic regime may be neutrally stable, it represents the universal small- 
scale asymptotics of steady forced turbulence. Experiments (Rutgers, 1998; Paret et ai, 
1999; Jullien et ai, 2000) and numerical simulations (Borue, 1993; Gotoh, 1998; Bowman 
et ai, 1999) are compatible with that conclusion. 

The physics of the enstrophy cascade is thus basically the same as that for a passive 
scalar: a fluid blob embedded into a larger-scale velocity shear is extended along the direc- 
tion of a positive strain and compressed along its negative eigendirection; such stretching 
provides for the vorticity flux toward the small scales, with the rate of transfer proportional 
to the strain. The vorticity rotates the blob decelerating the stretching due to the rotation 
of the axes of positive and negative strain. One can show that the vorticity correlators are 
indeed solely determined by the influence of larger scales (that give exponential separation 
of the fluid particles) rather than smaller scales (that would lead to a diffusive growth as 
the square root of time). The subtle differences from the passive scalar case come from the 
active nature of the vorticity. Consider for example the relation (135) expressing the fact 
that the correlation function of a passive scalar is essentially the time spent by the particles 
at distances smaller than L. The passive nature of the scalar makes Lagrangian trajecto- 
ries independent of scalar pumping which is crucial in deriving (135) by two independent 
averages. For an active scalar, the two averages are coupled since the forcing affects the 
velocity and thus the Lagrangian trajectories. In particular, the statistics of the forcing 
along the Lagrangian trajectories <p (R(t)) involves nonvanishing multipoint correlations at 
different times. Falkovich and Lebedev (1994) argued that, as far as the dominant logarith- 
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mic scaling of the correlation functions is concerned, the active nature of vorticity simply 
amounts to the following: the field can be treated as a passive scalar, but the strain acting 
on it must be renormalized with the scale. Their arguments are based on the analysis of 
the infinite system of equations for the variational derivatives of the vorticity correlation 
functions with respect to the pumping and the relations between the strain and the vorticity 
correlation functions. The law of renormalization is then established as follows, along the 
line suggested earlier by Kraichnan (1967, 1971, 1975). From (12), one has the dimensional 
relation that times behave as u; _1 ln(L/r). Furthermore, by using the relation (135) for 
the vorticity correlation function, one has {ujuj) oc P 2 x time. Combining the two previous 
relations, the scaling 00 ~ [P 2 ln^/r)] 1 / 3 follows. The consequences are that the distance 
between two fluid particles satisfies: \n(R/r) ~ P^ 2 t 3 / 2 , and that the pair correlation func- 
tion {U1U2) ~ (P2 hi(L/r 12 )) 2 / 3 . Note that the fluxes of higher powers tu 2n are not constant 
in the inertial range due to the same phenomenon of "distributed pumping" discussed in 
(vi) of Sect. III.C.l. The vorticity statistics is thus determined by the enstrophy production 
rate alone. 

It is worth stressing that the logarithmic regime described above is a small-scale asymp- 
totics of a steady turbulence. Depending on the conditions of excitation and dissipation, 
different other regimes can be observed either during an intermediate time or in an interme- 
diate interval of scales. First, a constant friction that provides for the velocity decay rate a 
prevails (if present) over viscosity at scales larger than (u/a) 1 ^ 2 . At such scales, the vortic- 
ity correlation functions are expected to behave as power laws rather than logarithmically, 
very much like for the passive scalar as described in (Chertkov, 1999) and in Section III.E. 
Indeed, the advective and the friction terms v • Vw and —av have again the same dimension 
for a smooth velocity. Nontrivial scaling is therefore expected, including for the second-order 
correlation function (and hence for the energy spectrum). The difficulty is of course that the 
system is now nonlinear and exact closed equations, such as those in Section III.E, are not 
available. For theoretical attempts to circumvent that problem by some approximations, 
not quite controlled yet, see (Bernard, 2000; Nam et al, 2000). Second, strong large-scale 
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vortices often exist with their (steeper) spectrum masking the enstrophy cascade in some 
intermediate interval of scales (Legras et al, 1988). 



2. On the energy cascades in incompressible turbulence 

The phenomenology of the energy cascade suggests that the energy flux e is a major 
quantity characterizing the velocity statistics. It is interesting to understand the difference 
between the direct and the inverse energy cascades from the Lagrangian perspective. The 
mean Lagrangian time derivative of the squared velocity difference is as follows 

The right hand side coincides with minus twice the flux and this gives the Lagrangian inter- 
pretation of the flux relations. In the 2d inverse energy cascade, there is no energy dissipative 
anomaly and the right hand side in the inertial range is determined by the injection term 
4(/ • v). The energy flux is negative (directed upscale) and the mean Lagrangian derivative 
is positive. On the contrary, in the 3d direct cascade the injection terms cancel and the 
right hand side becomes equal to — Av{(Vv) 2 ). The mean Lagrangian derivative is nega- 
tive while the flux is positive (directed downscale). This is natural as a small-scale stirring 
causes opposite effects with respect to a small-scale viscous dissipation. The negative sign 
of the mean Lagrangian time derivative in 3c? does not contradict the fact that any couple 
of Lagrangian trajectories eventually separates and their velocity difference increases. It 
tells however that the squared velocity difference between two trajectories generally behaves 
in a nonmonotonic way: the transverse contraction of a fluid element makes initially the 
difference between the two velocities decrease, while eventually the stretching along the 
trajectories takes over (Pumir and Shraiman, 2000). 

The Eulerian form of (207) is the generalization of 4/5-law (2) for the <i-dimensional 
case: ((A r v) 3 ) = —12er/d(d + 2) if A r v is the longitudinal velocity increment and e is 
positive for the direct cascade and negative for the inverse one. Since the average velocity 
difference vanishes, a negative ((A r v ) 3 ) means that small longitudinal velocity differences 
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are predominantly positive, while large ones are negative. In other words, in 3d if the 
longitudinal velocities of two particles differ strongly then the particles are likely to attract 
each other; if the velocities are close, then the particles preferentially repel each other. 
The opposite behavior takes place in 2d, where the third-order moment of the longitudinal 
velocity difference is positive. Another Lagrangian meaning of the flux law in 3c? can be 
appreciated by extrapolating it down to the viscous interval. Here, A r v ar and the 
positivity of the flux is likely to be related to the fact that the negative Lyapunov exponent 
is the largest one (in absolute value) in 3d incompressible turbulence. 

If one assumes (after Kolmogorov) that e is the only pumping-related quantity that de- 
termines the statistics then the separation between the particles R 12 = R(t; 7*1) — R(t; r 2 ) 
has to obey the already mentioned Richardson law: (R\ 2 ) ^£ 3 . The equation for the 
separation immediately follows from the Euler equation (205): d?Ri 2 = f 12 — VPi 2 . The 
corresponding forcing term f 12 = f(R(t; 7*1)) — f(R(t; r 2 )) has completely different prop- 
erties for an inverse energy cascade in 2d than for a direct energy cascade in 3d. For the 
former, i? 12 in the inertial range is much larger than the forcing correlation length. The 
forcing can therefore be considered short-correlated both in time and in space. Was the 
pressure term absent, one would get the separation growth: (R\ 2 )/et 3 = 4/3. The exper- 
imental data by Jullien et al. (1999) give a smaller numerical factor ~ 0.5, which is quite 
natural since the incompressibility constrains the motion. What is however important to 
note is that already the forcing term prescribes the law {R\ 2 ) oc t 3 consistent with the scaling 
of the energy cascade. Conversely, for the direct cascade the forcing is concentrated at the 
large scales and f\ 2 oc R\ 2 in the inertial range. The forcing term is thus negligible and 
even the scaling behavior comes entirely from the advective terms (the viscous term should 
be accounted as well). Another amazing aspect of the 2d inverse energy cascade can be 
inferred if one considers it from the viewpoint of vorticity. Enstrophy is transferred toward 
the small scales and its flux at the large scales (where the inverse energy cascade is taking 
place) vanishes. By analogy with the passive scalar behavior at the large scales discussed in 
Sect. III.C.l, one may expect the behavior (U1U2) oc r\ 2 a ~ d ) where a is the scaling exponent 
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of the velocity. The self-consistency of the argument dictated by the relation u = V x v 
requires 1 — a — d = 2a — 2 which indeed gives the Kolmogorov scaling a — 1/3 for d — 2. 
Experiments (Paret and Tabeling, 1997, 1998) as well as numerical simulations (Smith and 
Yakhot, 193; Boffetta et al, 2000) indicate that the inverse energy cascade has a normal Kol- 
mogorov scaling for all measured correlation functions. No consistent theory is available yet, 
but the previous arguments based on the enstrophy equipartition might give an interesting 
clue. To avoid misunderstanding, note that in considering the inverse cascades one ought 
to have some large-scale dissipation (like bottom and wall friction in the experiments with 
a fluid layer) to avoid the growth of condensate modes on the scale of the container. An- 
other example of inverse cascade is that of the magnetic vector potential in two-dimensional 
magnetohydrodynamics, where the numerical simulations also indicate that intermittency 
is suppressed (Biskamp and Bremer, 1994). The generality of the absence of intermittency 
for inverse cascades and its physical reasons is still an open problem. The only inverse 
cascade fully understood is that of the passive scalar in Section III. E, where the absence of 
anomalous scaling was related to the uniqueness of the trajectories in strongly compressible 
flow. That explanation applies neither to 2d Navier-Stokes nor to magnetohydrodynamics 
since the scalar is active in both cases. Qualitatively, it is likely that the scale-invariance 
of an inverse cascade is physically associated to the growth of the typical times with the 
scale. As the cascade proceeds, the fluctuations have indeed time to get smoothed out and 
not multiplicatively transferred as in the direct cascades, where the typical times decrease 
in the direction of the cascade. 

An interesting phenomenological Lagrangian model of 3d turbulence based on the con- 
sideration of four particles was introduced by Chertkov, Pumir et al, (1999). As far as an 
anomalous scaling observed in the 3d energy cascade is concerned, the primary target is to 
understand the nature of the statistical integrals of motion responsible for it. Note that 
the velocity exponent a 3 = 1 and experiments demonstrate that a p — > ap as p — > with a 
exceeding 1/3 beyond the measurement error, see (Sreenivasan and Antonia, 1997) and the 
references therein. The convexity of a p means then that o<i > 2/3. In other words, already 
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the pair correlation function should be determined by some nontrivial conservation law (like 
for magnetic fields in Sect. III. C. 3). 

V. CONCLUSIONS 

This review is intended to bring home to the reader two main points: the power of the 
Lagrangian approach to fluid turbulence and the importance of statistical integrals of motion 
for systems far from equilibrium. 

As it was shown in Sects. II and III, the Lagrangian approach allows for a systematic 
description of most important aspects of particle and field statistics. In a spatially smooth 
flow, Lagrangian chaos and exponentially separating trajectories are generally present. The 
associated statistics of passive scalar and vector fields is related to the statistics of large 
deviations of stretching and contraction rates in a way that is well understood. The the- 
ory is quite general and it finds a natural domain of application in the viscous range of 
turbulence. The most important open problem here seems to be the understanding of the 
back-reaction of the advected field on the velocity. That would include an account of the 
buoyancy force in inhomogeneously heated fluids, the saturation of the small-scale magnetic 
dynamo and the polymer drag reduction. In nonsmooth velocities, pertaining to the inertial 
interval of developed turbulence, the main Lagrangian phenomenon is the intrinsic stochas- 
ticity of the particle trajectories that accounts for the energy depletion at short distances. 
This phenomenon is fully captured in the Kraichnan model of nonsmooth time-decorrelated 
velocities. To exhibit it for more realistic nonsmooth velocities and to relate it to the hy- 
drodynamical evolution equations governing the velocity field remains an open problem. 
The spontaneous stochasticity of Lagrangian trajectories enhances the interaction between 
fluid particles leading to intricate multi-particle stochastic conservation laws. Here, there 
are open problems already in the framework of the Kraichnan model. First, there is the 
issue of whether one can build an operator product expansion, classifying the zero modes 
and revealing their possible underlying algebraic structure, both at large and small scales. 
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The second class of problems is related to a consistent description of high-order moments 
of scalar, vector and tensor fields. In the situations where the amplitudes of the fields are 
growing, this would be an important step towards a description of feedback effects. 

Our inability to derive the Lagrangian statistics directly from the Navier-Stokes equations 
of motion for the fluid particles is related to the fact that the particle coupling is strong 
and nonlocal due to pressure effects. Some small parameter for perturbative approaches, 
like those discussed for the Kraichnan model, has often been sought for. We would like 
to stress, however, that most strongly coupled systems, even if local, are not analytically 
solvable and that not all measurable quantities may be derived from first principles. In 
fluid turbulence, it seems more important to reach basic understanding of the underlying 
physical mechanisms, than it is to find out the numerical values of the scaling exponents. 
Such an understanding has been achieved in passive scalar and magnetic fields through the 
statistical conservation laws. We consider the notion of statistical integrals of motion to be 
of central importance for fluid turbulence and general enough to apply to other systems in 
non-equilibrium statistical physics. Indeed, non-dimensionally scaling correlation functions 
appearing in such systems should generally be dominated by terms that solve dynamical 
equations in the absence of forcing (zero modes). As we explained throughout the review, 
for passively advected fields, such terms describe conservation laws that are related to the 
geometry either of the configuration of particles (for the scalar) or of the particle-plus-field 
configurations (for the magnetic field). It is a major open problem to identify the appropriate 
configurations for active and nonlocal cases. New particle-tracking methods (La Porta et 
al, 2001) open promising experimental possibilities in this direction. An investigation of 
geometrical statistics of fluid turbulence by combined analytical, experimental and numerical 
methods aimed at identifying the underlying conservation laws is a challenge for future 
research. 
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FIGURES 

FIG. 1. An illustration of the breakdown of the Lagrangian flow in spatially non-smooth flows: 
infinitesimally close particles reach a finite separation in a finite time. The consequence is the cloud 
observed in the figure. The particles evolve in a fixed realization of the velocity field and in the 
absence of any molecular noise. 

FIG. 2. The implosion of Lagrangian trajectories in a strongly compressible flow. Particles that 
are initially released uniformly across a sizable span of the interval are compressed and tend to 
produce a singular density field. 

FIG. 3. An example of Lagrangian trajectories of three particles. The probability density of the 
positions R, conditional to the r's, is described by the PDF p(r, s;R,t\ v) (in a fixed realization 
of the velocity). Its average over the statistics of the velocity field gives the Green functions 
V(r; E; t - s). 

FIG. 4. The contour lines of a three particle zero mode as a function of the shape of the triangle 
defined by the particles. 

FIG. 5. The fourth-order anomalous exponent 2^— C4 of the scalar field rathe roughness parameter 
£ of the velocity field in the Kraichnan model. The circles and the stars refer to the three- 
dimensional and the two-dimensional cases, respectively. The dashed lines are the perturbative 
predictions for small £ and 2 — £ in 3d. 

FIG. 6. A typical snapshot of a scalar field transported by a turbulent flow. 

FIG. 7. The PDF's V(A r 0) of the scalar increments A r 9 = 0(r) - 9(0) for three values of r inside 
the inertial range of scales, multiplied by the factor r"^. The observed collapse of the curves 
implies the saturation of the scaling exponents of the scalar structure functions. 
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FIG. 8. Geometric construction of the Hopf-Cole inviscid solution of the ID Burgers equation. The 
inverted parabola C — ^(r — r') 2 is moved upwards until the first contact point with the profile 
of the initial potential </>(r',0). The corresponding height C gives the potential 4>(r,t) at time t. 
Shocks correspond to positions r where there are several contact points r', as for the first parabola 
on the left. 
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APPENDIX A: REGULARIZATION OF STOCHASTIC INTEGRALS 



For the 5-correlated strain, the equation (12) becomes a stochastic differential equation. 
Let us present here a few elementary facts about such equations for that simple case. The 
differential equation (12) is equivalent to the integral equation 

t 

R(t) = R(0) + f <t(s) ds R(s) , (Al) 
o 

where R = R± 2 . The right hand side involves a stochastic integral whose distinctive feature 
is that a(t)dt is of order (dt) 1 ^ 2 , as indicated by the relation ^ J * <7y(s) ds^j ^ oc t. Such 
integrals require second-order manipulations of differentials and, in general, are not unam- 
biguously defined without the choice of a defining convention. The most popular are the 
Ito, the Stratonovich and the anti-Ito ones. Physically, different choices reflect finer details 
of the strain correlations wiped out in the white- noise scaling limit, like the presence or the 
absence of time-reversibility of the velocity distribution. The Ito, Stratonovich and anti-Ito 
versions of the stochastic integral in (Al) are given by the limits over partitions of the time 
interval of different Riemann sums: 



ds R(s) 



tn+l 

lim / &( s ) ds R(t n ) , 
n t n 

*n+l 

lim E / *{s) ds l[R(t n ) + Jl(f n+ i)] , (A2) 

t n +l 

lim E / cr(s) ds R(t n+1 ) , 



respectively, where = t < t\ < ■ ■ ■ < tN — t. It is not difficult to compare the different 
choices. For example, the difference between the second and the first one is 

ilim^ / a(s) ds [R(t n+1 ) - R(t n )\ = \ lim ^ ( / a(s)ds) R(t n ) 

= JCR{s)ds, (A3) 
o 

where Ca = Cijje (sum over j) and the last equality is a consequence of the Central Limit 

Theorem that suppresses the fluctuations in lim J2 ( / a(s) ds) . Similarly, the difference 

n v t n ' 

between the anti-Ito and the Ito procedure is twice the latter expression. In other words, the 
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Stratonovich and anti-Ito versions of (12) are equivalent, respectively, to the Ito stochastic 
equations 

([a(t)dt + Cdt]R, 
dR = I (A4) 
{ [a(t) dt + 2C dt] R . 

Given the stochastic equation (12) with a fixed convention, its solution can be obtained 
by iteration from (Al) and has the form (13) with W(t) given by (16) and the integrals 

I t s , 

interpreted with the same convention. Note that the value of < / a(s) ds J a(s') ds') depends 

x o o 1 

on the choice of the convention: it vanishes for the Ito one and is equal to C t for the 
Stratonovich and twice that for the anti-Ito ones. The conventions are clearly related to the 
time-reversibility of the finite-correlated strain before the white-noise scaling limit is taken. 
For example, time-reversible strains have even 2-point correlation functions and produce the 
Stratonovich value for the above integrals in the white- noise scaling limit. Most of these 
stochastic subtleties may be forgotten for the incompressible strain where C = as follows 
from (30) so that the difference between different conventions for (12) disappears. 

We shall often have to consider functions of solutions of stochastic differential equations 
so one should be aware that the latter behave under such operation in a somewhat peculiar 
way. For the Ito convention, this is the content of the so-called Ito formula which results 
from straightforward second-order manipulations of the stochastic differentials and takes in 
the case of (12) the form: 

df(R) = (a(t)dtR)-Vf(R) + CijuR? R l WiW k f\R)dt. (A5) 

Note the extra second-order term absent in the normal rules of differential calculus. The 
latter are, however, preserved in the Stratonovich convention. Note that the latter difference 
may appear even for the incompressible strain. 

More general stochastic equations may be treated similarly. For example, (5) for Kraich- 

nan velocities may be rewritten as an integral equation involving the stochastic integral 

t 

J v(R(s), s) ds. The latter is defined as 
o 
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lim E / ^[v(R(t n ),s) + v(R(t n+1 ),s)]ds (A6) 

n , 

in the Stratonovich convention, with the last t n+ i (t n ) replaced by t n (t n+ i) in the (anti-)Ito 
one. The difference is 

-, tn + l tn + 1 

±\ lim 2 / / (u(fl(t n ),s') • V) v(il(* n ),s) rfs'rfs = ±V^(0)t, (A7) 

72 in in 

where the left hand side was replaced by its mean by virtue of the Central Limit Theorem. 
The last term vanishes if the velocity 2-point function is isotropic and parity invariant. 
Hence the choice of the convention is unimportant here even for the compressible velocities. 

It may seem bizarre that the choice of convention in compressible velocities does not 
matter for individual trajectories but it does for the equation (12) which describes the 
evolution of small trajectory differences. It is not difficult to explain this discrepancy (Horvai, 
2000). The stochastic equation for R 12 = R leads (in the absence of noise) to the integral 
equation 

t 

R(t) = R(0) + J [v(R(s) + R 2 (s),s) - v(R 2 (s),s)} ds . (A8) 
o 

The difference between the Stratonovich and the Ito conventions for the latter integral is 

-, tn+l tn+1 

5 Um 2 / / {[v(R(t n ) + R 2 (t n ),s')-v(R 2 (t n ),s')]-V}v(R(t n ) + R 2 (t n ),s) 

fl tn tn 
t n +l t n + l 

+\ lim I I v{R 2 {t n ),s') ■ V)[v(R(t n ) + R 2 (t n ), s) - v(R 2 (t n ),s)} ds 

n in tn 

t t 
= J[V^{R{s)) - V^'(O)] ds + /[V^'(O) - V 3 D^{R{s))] ds = , (A9) 
o o 

where the two first lines that cancel each other are due to the time dependence of R and 
R 2 , respectively. If we replace R by eR then, when e — > 0, the right hand side of (A8) 
is replaced by the right hand side of (Al) if we use the Ito convention for the stochastic 
integrals. The similar limiting procedure applied to the (vanishing) difference (A9) does not 
reproduce the difference between the values of the integral (Al) for different conventions. 
The latter corresponds to the limit of the first line of (A9) only and does not reproduce the 
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limit of the term of (A9) due to the independence. The approximation (12) is then valid 
only within the Ito convention. 

As for the PDF (43), it may be viewed as a (generalized) function of the trajectory 
R(t). The advection equation (45) results then from the equation for R(t) by applying the 
standard rules of differential calculus which hold when the Stratonovich convention is used. 
The Ito rules produce the equivalent Ito form of the equation with an additional second 
order term containing the eddy diffusion generator D V 2 R . 
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